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ABSTRACT 

An  orbiting  detector  of  infrared  (IR)  energy  may  be  used  to  detect  the  rocket 
plumes  generated  by  ballistic  missiles  during  the  powered  segment  of  their  trajectory  By 
measuring  angular  directions  of  the  detections  over  several  observations,  the  trajectory 
properties,  launch  location,  and  impact  area  may  be  estimated  using  a  nonlinear  least- 
squares  iteration  procedure.  Observations  from  two  or  more  sensors  may  be  combined  to 
form  stereoscopic  lines-of-sight  (LOS),  increasing  the  accuracy  of  the  estimation 
algorithm  The  focus  of  this  research  has  been  to  develop  a  computer-model  of  an 
estimation  algorithm,  and  determine  what  parameter,  or  combination  of  parameters  will 
significantly  affect  on  the  error  of  the  tactical  parameter  estimation  This  model  is  coded 
in  MATLAB,  and  generates  observation  data,  and  produces  an  estimate  for  time,  position, 
and  heading  at  launch,  at  burnout,  and  calculates  an  impact  time  and  position.  The  effects 
of  time  errors,  LOS  measurement  errors,  and  satellite  position  errors  upon  the  estimation 
accuracy  were  then  determined  using  analytical  and  Monte  Carlo  simulation  techniques. 
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A  matrix  of  partial  derivatives  of  focal  plane  coordinates  with  respect  to 
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I.  INTRODUCTION 

The  TRW-built  Defense  Support  Program  (DSP)  satellites  have  been  the 
spaceborne  segment  of  NORADs  Tactical  Warning  and  Attack  Assessment  system  since 
the  early  1970's  Using  infrared  detectors  that  sense  the  heat  from  missile  plumes  against 
the  Earth  background,  these  orbiting  sentries  detect  ballistic  missile  launches  DSP  also 
detects  nuclear  detonations  in  support  of  Nuclear  Test  Ban  monitoring  The  DSP  system 
provides  near  real-time  detection  information  in  support  of  DOD's  tactical  warning  and 
attack  assessment  mission  and  is  supported  by  a  network  of  fixed  and  mobile  ground 
stations  that  process  and  disseminate  information  to  military  commanders  worldwide  The 
Cold  War  mission  of  the  DSP  system  was  to  detect  massive  intercontinental  ballistic 
missile  (ICBM)  attacks  United  States'  response  to  such  an  attack  only  required  timely 
and  unambiguous  warning  Missile  flight  times  were  much  longer  than  the  time  required 
to  launch  a  retaliatory  attack.  Precise  radar  tracks  could  be  established  with  enough  time 
to  mitigate  effects  as  much  as  technology  allowed. 

The  current  combat  environment  demands  much  more  from  launch  detection 
satellites  The  present  threat  is  from  tactical  ballistic  missiles  (TBM's),  which  exhibit 
much  cooler  burn  and  shorter  thrust  times,  and  possibly  more  depressed  trajectories  than 
those  exhibited  by  ICBMs  TBM's  can  be  launched  from  almost  anywhere  within  a  large 
geographical  area  of  interest,  with  lofted  or  depressed  trajectories  There  may  be  no  other 
sensors  to  quantify  impact  parameters  in  time  to  employ  anti-ballistic  missile  (ABM) 
weapons  or  alert  potential  victims  within  the  impact  zone. 

For  budgetary  reasons,  the  United  States  is  forced  to  use  existing  launch  detection 
satellites  to  counter  the  TBM  threat  into  the  beginning  of  the  next  century.  [Ref  1]  More 
rapid  extraction  of  more  precise  information  from  these  existing,  technology-limited 
satellites  must  be  accomplished  in  the  interim  of  acquiring  the  next-generation  TBM 
detection  satellite  system  Anti-ballistic  missile  (ABM)  systems,  such  as  Patriot  or  Aegis, 
may  be  employed  as  ABM  umbrellas  in  a  tactical  area  of  interest,  provided  they  receive 
precise  and  timely  cueing  from  these  detection  systems. 


This  raises  the  question:  "How  accurate  is  the  information  provided  by  DSP9"  To 
begin  to  answer  this  question,  an  engineering  error  analysis  must  be  performed  to 
determine: 

•  What  are  the  errors  present  in  the  detection  system? 

•  Which  errors  have  the  greatest  effect  upon  the  accuracy  of  DSP  output 
information9 

•  What  are  the  effects  of  the  errors  on  the  results  (individually  and  collectively)? 

By  modeling  the  algorithm  used  by  the  tactical  warning  system  to  determine  trajectory 
elements  and  predict  impact  zone  location  of  a  detected  TBM,  it  is  then  possible  to 
introduce  the  inherent  errors  separately  into  the  model,  and  then  study  their  effects  upon 
the  results.  The  relative  magnitudes  of  their  effects  upon  the  output  and  their  effects  upon 
the  results  are  then  determined.  This  is  done  both  qualitatively  and  numerically 
(statistically),  and  the  results  are  compared. 


II.  BACKGROUND 


The  DSP  system  consists  of  one  or  more  satellites  in  geosynchronous  orbit  and 
one  or  more  ground  receiving  stations  A  DSP  satellite  consists  of  a  bus  (spacecraft) 
segment  and  a  payload  (sensor)  segment.  The  satellite  is  10  meters  long,  7  meters  in 
diameter,  and  weighs  over  2300  kilograms.  [Ref.  2] 
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Figure  1    Diagram  of  DSP  Satellite  [From  Ref.  2] 

The  spacecraft  segment  receives  and  transmits  all  commands  and  data,  provides 
attitude,  spin-rate,  and  station-keeping  functions,  and  provides  and  distributes  satellite 
power.  Solar  cells  mounted  on  the  spacecraft's  cylindrical  body  and  on  four  solar  panels 
mounted  opposite  the  sensor  generate  the  electrical  power.  The  sensor  collects  IR  energy, 
provides  onboard  (initial)  data  processing,  and  supplies  precise  orientation  data 


The  satellite  is  placed  in  a  circular,  equatorial,  geosynchronous  orbit  (Titan/IUS), 
oriented  so  that  the  telescope  is  pointed  toward  the  earth,  and  spun  along  its  longitudinal 
axis  at  six  rpm.  This  configuration  is  called  a  "yaw  spinner".  The  axis  of  rotation  of  the 
satellite  is  normal  to  the  earth's  surface  (nadir-pointing).  The  satellite's  spin  allows  the 
sensor  to  regularly  scan  the  earth  and  distributes  the  thermal  load  uniformly.  The  spin  also 
allows  one  set  of  attitude  control  components  to  effect  two-axis  earth  pointing  control 
(roll  and  pitch)  by  time-sharing  A  counter-rotating  momentum  wheel  keeps  the  net  spin 
momentum  near  zero,  thereby  minimizing  the  gyroscopic  effect  due  to  coupling  between 
the  spin  motion  and  the  orbit  rate  with  minimum  fuel  expenditure 

The  IR  telescope  is  tilted  from  the  spin  axis,  so  that  the  photo-electric  cell  (PEC) 
array  covers  the  radius  of  the  earth.  As  the  satellite  rotates,  the  entire  surface  of  the  earth 
within  the  FOV  is  scanned  by  the  IR  detector,  shown  in  Figures  2  and  3 . 
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Figure  2    Photoelectric  Cell  Array  [From  Ref  2] 
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Figure  3    Focal  Plane  Array  Scanning  FOV  [From  Ref  3] 


Detection  of  IR  sources  is  accomplished  with  the  telescope  and  PEC-array 
portions  of  the  sensor.  The  signal  electronics  onboard  the  satellite  partially  discriminate 
the  targets  from  the  background  and  the  ground  station  completes  the  task.  Location  of 
the  IR  sources  is  derived  from  orientation  of  the  IR  telescope  line-of-sight  relative  to  the 
earth's  surface  Star  sensors  augment  the  sensor  data  for  precise  determination  of  sensor 
pointing  direction. 

The  PEC-array  of  the  IR  detector  is  mounted  with  the  nadir  end  at  the  center  of 
rotation  of  the  telescope  (see  Figure  3).  This  array  contains  over  6000  detector  cells  that 
are  sensitive  to  energy  in  the  infrared  wavelengths.  As  the  PEC  array  scans  the  FOV,  a 
cell  passing  across  an  IR  source  will  generate  a  voltage  with  an  amplitude  proportional  to 
the  signal  intensity.  This  voltage  signal  is  termed  an  IR  return,  and  is  transmitted  to 
ground  processing  stations  after  amplification  and  background  filtering.  A  simplified 
diagram  of  the  data  distribution  network  is  shown  in  Figure  4 


GROUND  DATA 


Figure  4    Data  Distribution  Diagram  [Ref  4] 
DSP  has  ground  stations  globally  positioned  to  receive,  process,  and  report 
mission  data  to  users.     The  Satellite  Control  Facility  tracks  the  satellite,  transmits 
commands,  receives  mission  and  satellite  status  data,  and  forwards  the  data  to  the  Data 


Reduction  Center  (DRC)  for  further  processing  The  DRC  processes  the  data,  extracts 
significant  returns  associated  with  missile  launches  and  nuclear  detonations,  and  generates 
event  messages  for  transmission  to  the  users  over  the  Ground  Communications  Network 
[Ref.  2] 


III.  THE  ALGORITHM 

{ Thomas  Jerardi  has  graciously  given  his  permission  to  adapt  and  modify  his 
unpublished  paper,  "TALON  SHIELD  ALERT  State  Vector  Estimation",  which  is  the 
basis  for  the  majority  of  this  chapter. }  [Ref.  5] 

The  detection  of  a  TBM  launch  is  the  starting  point  for  any  ballistic  missile  de- 
fense. DSP  does  this  by  detecting  the  IR  radiation  emitted  by  the  exhaust  plume  of  a 
launching  missile.  With  detections  by  two  or  more  spacecraft  (stereo  observations),  trian- 
gulation  of  lines  of  sight  can  be  used  to  more  accurately  estimate  the  boost-phase  trajec- 
tory, which  is  then  used  to  calculate  launch  position,  state  vector  (position  and  velocity)  at 
missile  engine  burnout,  and  impact  position  Real-time  knowledge  of  the  launch  position 
allows  targeting  of  the  launcher.  Cueing  for  ABM  weapons  systems,  such  as  Patriot  or 
Aegis,  requires  timely  and  accurate  trajectory  information,  which  can  be  propagated  from 
knowledge  of  the  state  vector  at  burnout  Impact  time  and  position  data  is  extrapolated 
from  the  state  vector  at  burnout,  and  may  be  used  for  warning  personnel  within  the  target 
area.  Therefore,  the  quality  (accuracy)  of  the  trajectory  estimation  process  is  of  para- 
mount importance  to  Tactical  Ballistic  Missile  Defense  (TBMD)  Understanding  the 
algorithms  and  equations  employed  in  the  estimation  process  is  necessary  to  assess  the 
quality  of  the  estimated  parameters. 

The  tactical  parameter  estimation  process  is  composed  of  several  tasks:  initial 
estimate  of  the  tactical  parameters,  nonlinear  least-squares  estimation  (refinement)  of 
tactical  parameters,  burn-out  time  estimation,  state  vector  generation,  and  impact  point 
calculation    The  tactical  parameters  are: 

•  To  =  Time  of  launch 

•  L  -Loft 

•  <J>o  =  Launch  point  geodetic  latitude 

•  Xo  =  Launch  point  geodetic  longitude 

•  ho  =  Launch  point  height  above  WGS-84  ellipsoid 

•  do  =  Flight  trajectory  azimuth  (true  heading) 

These  quantities  will  be  explained  in  more  detail  in  the  following  sections 


Observational  data  are  used  to  calculate  initial  estimates  of  tactical  parameters  and 
to  choose  the  appropriate  TBM  profile  from  a  database.  The  profile  trajectory  is  then 
used  to  calculate  expected  observations,  which  are  then  compared  to  the  actual  observa- 
tions. The  differences  are  used  to  calculate  changes  to  the  initial  estimates,  and  a  least- 
squares  iteration  is  run  until  the  differences  between  expected  and  actual  observations  are 
sufficiently  small.  The  result  is  an  accurate  determination  of  the  TBM's  launch  and  trajec- 
tory parameters.  An  estimate  of  burnout  time  is  then  made,  and  the  calculated  trajectory 
is  extrapolated  to  generate  the  state  vector  at  burnout.  This  allows  the  computation  of 
impact  time  and  position  using  simple  ballistic  trajectory  equations  Each  of  these  tasks  is 
described  in  detail  in  the  following  sections 

A.  OBSERVATIONAL  DATA 

The  starting  point  is  the  observational  data,  which  are  measurements  of  IR  radiant 
intensities  (IR  returns)  as  a  function  of  time  from  each  of  the  approximately  6000  focal 
plane  elements  in  the  DSP  satellite's  sensor  Attitude  information  (star  sensor  measure- 
ments, jet-firing  data,  momentum  wheel  tachometer  data)  is  also  included  in  the  telemetry 
stream  and  used  to  determine  spacecraft  attitude  history.  The  satellite's  position  is  de- 
termined by  the  Air  Force  Satellite  Control  Network  (AFSCN).  A  global  perspective  of 
the  observation  geometry  is  shown  in  Figure  5. 


Figure  5.  Observation  Geometry  [From  Ref.  5] 

10 


The  discrimination  of  TBM  IR  returns  from  earth  background  is  a  complex  proc- 
ess It  will  be  assumed  to  have  been  effectively  done,  and  that  a  table  of  observations 
corresponding  to  a  single,  boosting  TBM  is  available    An  example  is  shown  in  Table  1 


Index 

Time 

(sec) 

S/C 
ID 

Intens. 

Azimuth 
(rad) 

Elevation 
(rad) 

G.H.A. 

(rad) 

Dec. 
(rad) 

Radius 
(km) 

k 

T 

S/C 

I 

P 

1 

gha 

6 

R 

1 

129.36 

1 

14.0 

4061854 

0.115847 

0.174532 

0 

42164  17 

2 

130.30 

2 

23.0 

2427020 

0.094741 

1.221730 

0 

42164  17 

3 

135.44 

3 

32.0 

2.062181 

0.142310 

1.832595 

0 

42164.17 

4 

139.36 

1 

29.0 

4.062497 

0  115971 

0.174532 

0 

42164.17 

5 

140.30 

2 

29.0 

2.427264 

0.094753 

1.221730 

0 

42164.17 

6 

14544 

3 

40.0 

2.061969 

0.142405 

1.832595 

0 

42164  17 

7 

149.36 

1 

49.0 

4.063552 

0.116147 

0.174532 

0 

42164.17 

8 

150.30 

2 

60.0 

2.427660 

0.094746 

1.221730 

0 

42164.17 

9 

155.44 

3 

65.0 

2.061752 

0.142493 

1.832595 

0 

42164  17 

Table  1  Example  Observational  Data 
The  index,  k,  runs  from  1  to  n,  the  total  number  of  observations  (typically  less  than 
20),  and  is  used  as  a  subscript  for  the  remaining  symbols  to  denote  a  particular  observa- 
tion. Time,  Tk,  is  the  time  of  observation  measured  from  midnight  of  the  day  of  the  obser- 
vation, and  runs  between  0  and  86400  seconds.  Spacecraft  Identification,  S/C  ID,  is 
simply  a  notation  used  to  identify  which  satellite  is  making  which  observation  The  in- 
tensity, Ik,  is  the  radiant  intensity  of  the  IR  return,  and  is  used  to  select  the  type  of  TBM 
being  detected.  Azimuth  angle,  p\,  is  the  azimuth  of  the  line  of  sight  of  the  IR  return 
measured  clockwise  from  true  south,  and  has  a  range  of  0  to  2n  radians  (0°  to  360°). 
Elevation  angle,  r|k,  is  the  elevation  of  the  return  measured  from  nadir,  and  has  values 
between  0  and  0.175  radians  (0°  to  10°).  Satellite  position  is  given  in  spherical  coordi- 
nates. Greenwich  hour  angle,  ghak,  is  the  angle  between  the  satellite's  nadir  point  and  the 
prime  meridian,  measured  positive  east.  Declination  angle,  8k,  is  the  angle  above  or  below 
the  equator,  measured  positive  north  (Greenwich  hour  angle  and  declination  have  been 
adopted  instead  of  sub-satellite  point  longitude  and  latitude  to  avoid  confusion  with  TBM 
latitude  and  longitude.)  Radius,  Rk,  is  the  distance  from  the  satellite  to  the  earth's  center, 
measured  in  kilometers  Geosynchronous  satellites  typically  stay  within  a  few  degrees  of 
the  equator,  and  have  a  radius  of  approximately  42,164. 17  km. 
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B.  TBM  PROFILE 


A  TBM  profile  is  a  description  of  the  nominal  powered  flight  trajectory  of  the 
given  TBM    An  example  TBM  profile  is  given  in  Table  2. 


Time 

(sec) 

Intensity 

Altitude 
(km) 

Range 
(km) 

Time 

(sec) 

Intensity 

Altitude 
(km) 

Range 
(km) 

0 

36.0 

0.000 

0.000 

32 

60.6 

7.023 

3.195 

1 

36.3 

0.006 

0.000 

33 

62.4 

7469 

3.491 

2 

36.6 

0.026 

0.000 

34 

64.2 

7.928 

3.803 

3 

36.9 

0.058 

0.000 

35 

66.0 

8.402 

4.132 

4 

37.2 

0.103 

0.000 

36 

68.4 

8.890 

4.479 

5 

37.5 

0.163 

0.001 

37 

70.8 

9.393 

4.844 

6 

37.8 

0.235 

0.004 

38 

73.2 

9.911 

5.229 

7 

38.1 

0.322 

0.010 

39 

75.6 

10.444 

5.633 

8 

38.4 

0.423 

0.020 

40 

78.0 

10.992 

6.057 

9 

38.7 

0.537 

0.036 

41 

81.2 

11.556 

6.502 

10 

39.0 

0.666 

0.058 

42 

84.4 

12.136 

6.969 

11 

39.5 

0.809 

0.087 

43 

87.6 

12.732 

7.459 

12 

40.0 

0.965 

0.124 

44 

90.8 

13.345 

7.973 

13 

40.5 

1.136 

0.171 

45 

94.0 

13.975 

8.511 

14 

41.0 

1.321 

0.226 

46 

96.0 

14.622 

9.075 

15 

41.5 

1.520 

0.292 

47 

98.0 

15.288 

9.665 

16 

42.0 

1.733 

0.367 

48 

100.0 

15.972 

10.282 

17 

42.5 

1.962 

0.453 

49 

102.0 

16.675 

10.928 

18 

43.0 

2.204 

0.550 

50 

104.0  — 

17.397 

11.604 

19 

43.5 

2.460 

0.658 

51 

104.6 

18.140 

12.309 

20 

44.0 

2.731 

0.777 

52 

105.2 

18.904 

13.045 

21 

45.0 

3.015 

0.908 

53 

105.8 

19.690 

13.813 

22 

46.0 

3.312 

1.050 

54 

106.4 

20.499 

14.613 

23 

47.0 

3.623 

1.205 

55 

107.0 

21.332 

15.446 

24 

48.0 

3.948 

1.372 

56 

106.4 

22.190 

16.314 

25 

49.0 

4.286 

1.551 

57 

105.8 

23.075 

1.217 

26 

50.6 

4.637 

1.744 

58 

105.2 

23.986 

18.155 

27 

52.2 

5.001 

1.950 

59 

104.6 

24.925 

19.131 

28 

53.8 

5.378 

2.170 

60 

104.0 

25.894 

20.145 

29 

55.4 

5.769 

2.404 

61 

98.0 

26.894 

21.199 

30 

57.0 

6.174 

2.652 

62 

80.0 

27.925 

22.293 

31 

58.8 

6.591 

2.916 

62.5 

20.0 

28.450 

22.850 

Table  2.  Sample  TBM  Profile  [Ref  5] 


A  profile  consists  oflR  intensity  as  a  function  of  time,  nominal  (maximum  range 
trajectory)  vertical  and  horizontal  ranges  from  the  launch  point  as  functions  of  time,  and 
maximum  burn  time,  tmax  (62  5  seconds  in  the  example)  The  detected  radiant  intensity 
over  time  is  compared  to  TBM  profiles  in  a  data  base,  and  the  best  match  is  selected  as  the 
type  of  TBM  being  observed  This  selection  process,  called  "typing",  is  complex,  and  is 
also  assumed  to  have  been  done  correctly  A  more  computationally  efficient  (but  equiva- 
lent) representation  of  the  nominal  downrange  distance  and  altitude  profiles  can  be  found 
by  fitting  quartic  polynomials  to  the  data  in  the  profile  tables.  Use  of  the  polynomial  repre- 
sentation during  calculations  precludes  table  look-up  and  interpolation  problems  for  non- 
integer  times-of-flight    The  profile  distance  polynomial  has  the  form: 

dp  =a0  +a,t  +  a2t2  +a3t3  +a4t4  (1) 

and  the  profile  height  polynomial  is: 

hp=b0+b,t  +  b2t2+b3t3+b4t4  (2) 

where  t  is  time  of  flight,  and  ai  4  and  bi  4  are  coefficients  of  the  fourth-order  polynomial. 
These  coefficients  are  not  computed  within  this  algorithm,  but  are  assumed  to  be  known, 
exact,  correct,  and  available  from  the  typing  process,  i.e.,  perfect  a  priori  knowledge  of 
the  particular  observed  TBM's  nominal  trajectory.  Since  real-life  TBM's  do  not  fly  the 
profile  exactly,  using  the  profile  to  determine  the  trajectory  introduces  an  error  into  the 
algorithm  The  inherent  error  of  inexact  profile  information  is  ignored  in  this  analysis. 
Assuming  the  profile  to  be  exact  does  not  invalidate  the  error  analysis,  however. 

A  "loft"  parameter,  L,  is  used  to  account  for  trajectories  above  or  below  the  nomi- 
nal profile.  This  very  simple  model  for  a  ballistic  missile  trajectory  adjusts  the  profile  to 
give  actual  range: 

d  =  (l-1.5L)dp  (3) 

and  altitude: 

h  =  (l  +  L)hp  (4) 

Loft  varies  over  approximately  ±0.25.  Figure  6  is  a  plot  of  nominal  (L  =  0),  lofted  (L  = 
+0.25)  and  depressed  (L  =  -0.25)  trajectories. 


Lofted 


Nominal 


Depressed 


Range  (km) 

Figure  6.  Trajectory  Profiles  for  Three  Loft  Values 

C.  LINE  OF  SIGHT  PROJECTION 

The  least-squares  iteration  process  begins  with  an  initial  estimate  of  the  tactical 
parameters  In  order  to  make  an  estimate  for  the  initial  TBM  position  (A«.  <J>oX  the  focal 
plane  measurements  must  be  changed  into  a  LOS  between  the  satellites  and  the  TBM, 
which  points  to  the  position  of  the  TBM  on  the  globe  The  first  step  is  to  transform  the 
satellite  positions  (gha,  8,  R)  into  earth-centered  rotating  coordinates  (XYZ),  in  which  the 
X-axis  extends  through  the  prime  meridian  and  the  Z  axis  is  aligned  with  earth's  spin  axis, 
depicted  in  Figure  7. 
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Figure  7.  Coordinate  Frame  Illustration 
The  satellite  position  transformation  from  (gha,  5,  R)  to  (XYZ)  is: 
(Rk  cosCSJcostgha^ 


zsky 


Rk  cos(8k)sin(ghak) 
Rksin(5k) 


(5) 


where,  the  subscript,  s,  denotes  satellite  position,  and  the  sub-subscript  refers  to  the  k* 
observation  The  focal  plane  coordinates,  rjk  and  |3k,  are  first  transformed  to  the  satellite's 
Up-East-North  reference  frame,  (UEN),  and  then  must  also  be  transformed  into  (XYZ) 
The  focal  plane  coordinates,  can  be  visualized  with  Figure  8. 


Figure  8.  Focal  Plane  Coordinates 
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The  transformation  of  polar  (r)k,  (5k)  coordinates  into  Cartesian  (UEN)  coordinates 


-cos(rik) 
-sin(pk)sin(rik)  (6) 

-cos(pk)sin(rik) 

This  is  the  unit  line-of-sight  (LOS)  direction  vector  in  (UEN)  coordinates.  Next  transform 
into  (XYZ)  (refer  to  Figure  7): 

cos(ghak)    -sin(ghak)    0Tcos(8k)    0    -sin(8k) 
ek  =    sin(ghak)      cos(ghak)     0         0  1  0        ek  (7) 

0  0  l_[sin(5k)     0     cos(8k) 

For  the  initial  position  estimate,  it  is  sufficient  to  assume  that  the  line  of  sight 
terminates  on  the  surface  of  a  spherical  earth  with  an  effective  radius  of  reff  =  6371  km. 
The  Cartesian  coordinates  of  the  TBM  are: 


>\ 


hPkek 


(8) 


where  the  subscript,  T,  denotes  TBM  position,  sub-subscript  k  refers  to  the  k*  observa- 
tion, and  p  is  the  length  of  the  LOS  vector  as  shown  in  Figure  7.  The  Law  of  Sines  is  one 
of  several  methods  that  may  be  used  to  calculate  p.  The  geometry  of  the  problem  is 
shown  in  Figure  9. 


DSP 


Figure  9.  Initial  Position  Estimate  Geometry 
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The  angle  opposite  p  is 


Rk|sin(r|k) 


Xk  =  K~r\k  -sin 

With  some  trigonometric  identities  and  algebra,  the  expression  for  p  becomes 

'|Rk|sin(Tik)^ 


n     -r^sin(Xk) 
Kk    ~  •      ,         v 

sin(r|k) 
The  longitude  of  this  point  is: 


refr  sm 


rik  +  sin 


rcff 


sin(rik) 


and  the  latitude  is: 


A,k  =  tan 


(|>k  =  tan  ' 


-fe 


Vxik2  +yTk2 


(9) 


(10) 


(11) 


(12) 


D.  INITIAL  ESTIMATES 

The  core  of  the  estimation  process  is  a  nonlinear  least-squares  iterative  calculation 
of  the  tactical  parameters.  This  method  is  based  on  a  one-term  (linear)  Taylor  expansion 
of  the  focal  plane  observations  in  terms  of  the  six  tactical  parameters  This  procedure 
requires  an  initial  estimate  of  the  tactical  parameters  The  accuracy  of  the  initial  estimate 
does  not  effect  the  accuracy  of  the  final  estimate,  but  can  effect  the  convergence  time  of 
the  least-squares  process. 

The  initial  estimate  of  the  launch  position  and  heading  is  based  on  the  projected 
positions  (fa  and  Xk)  obtained  from  the  LOS  observations  The  lines-of-sight  from  each 
satellite  usually  do  not  intersect  at  a  single  point  (even  if  they  were  simultaneous  observa- 
tions), as  shown  in  Figure  10  The  initial  latitude  and  longitude  guesses  are  found  by 
calculating  the  average  position  from  the  first  LOS  observation  from  each  satellite  The 
assumption  is  made  that  if  one  satellite  can  see  the  TBM,  then  all  can  observe  it  This 
artificiality  is  not  true  in  the  physical  world 
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□  DSP  1 
A  DSP  2 
O       DSP  3 


O     □ 


46.6  46  7  46.8  46.9  47  0  47   1  47.2  47  3 

Longitude 

Figure  10    LOS  Projection  Scatter 
In  this  example,  three  satellites  observe  the  TBM,  so  the  first  LOS  projected  lati 
tudes  and  longitudes  (from  each  satellite)  are  averaged: 

.     (0)    _    <t>l   +<t>2   +<t>3 


X}  +X2  +x~ 


(13) 


(14) 


The  superscript, (0>,  denotes  the  initial  estimate    The  last  LOS  projected  positions  are  also 
averaged  to  obtain  the  approximate  final  position  observed 


*.-=- 


-Xa_l+Xm 


(15) 
(16) 


Equations  (13)  to  (16)  assume  that  the  first  and  last  three  observations  are  from  three 
different  satellites  If  only  two  satellites  observe  the  launch,  only  the  first  and  last  two 
observations  would  be  averaged. 

The  direction  of  the  last  position  from  the  first  is  the  initial  estimate  for  the  head- 

(0) 

ing,  ao  ': 

ar^-tan/'t^-C^cos^r),^,-^]  (17) 

The  function  tan2"'  is  the  two-argument  arctangent,  used  here  because  its  range  is  ±k. 
This  function  is  used  because  a  four-quadrant  answer  is  required 

The  first  estimate  for  T0  is  twenty  seconds  before  the  first  observation: 

To(0)  =  T,-20  (18) 

and  the  initial  guess  for  L  and  ho  are  both  zero: 

L<0>  =  0  (19) 

h0(0)   =  0  (20) 

E.  EXPECTED  POSITIONS  ALONG  THE  TRAJECTORY 

Given  the  six  (estimated)  tactical  parameters  at  launch,  the  expected  position  of 
the  TBM  along  its  boost  trajectory  can  be  determined  by  applying  the  TBM  profile  to  any 
observation  time,  and  in  particular.  Tk,  by  first  computing: 

tk=Tk-T0  (21) 

Equations  (1)  through  (4)  then  give  the  expected  altitude  (hk)  and  range  (dk)  from  the 
estimated  launch  point  when  t  =  tk: 

d,,k  =a0  +a,tk  +a2tk2  +a3tk3  +  a4tk4  (1) 

hpk  =b0+b,tk  +b2tk2  +b3tk3+b4tk4  (2) 

dk  =(1-1  5L)dP  (3) 

hk=(l  +  L)hPk  (4) 
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Next,  it  is  necessary  to  calculate  the  "earth-central  angle",  a  quantity  used  in  (23)  and 
(24): 


e4 


reff 


(22) 


0k  is  simply  the  angle  between  the  launch  point  and  the  position  of  the  TBM  with  the 
earth's  center  as  the  vertex    0  can  be  visualized  with  Figure  1 1 


ren        |     Launch 
position 


Figure  1 1 .  Earth-central  Angle 
Using  8k,  spherical  trigonometry  gives  the  equations  for  the  geodetic  latitude: 

cj)k    =    --cos-'[cos(0k)sin((|)o)  +  sin(0k)cos((t)o)cos(ao)]  (23) 


and  longitude: 


Xk   =  X0  +sin"' 


sin(9k)sin(a0) 


(24) 


cos(<|)k) 
The  altitude  is  simply  the  present  height  of  the  TBM  added  to  the  initial  height  at  launch: 

altk  =  h0+hk  (25) 

In  addition  to  the  geodetic  coordinates,  the  Cartesian  coordinates  are  also  required    Using 


re  =  6378.137  km,  and  f  = 


1 


298.257 
the  geocentric  latitude  of  the  TBM  is  calculated: 

4>'k  =  tan  '[(l-f)2tan((t)k)] 


,  the  WGS-84  values  for  the  flattening  of  the  earth, 


(26) 
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as  well  as  the  local  radius: 


(27) 


V(l-f)2cos2(fk)  +  sin2(<j>'k) 
These  values  are  then  transformed  into  earth-centered,  Cartesian  coordinates  of  the  TBM 
x -ik  =  [fbcu  «*C*\ ) +  ^k  C0S(<K  )]cos(X.k )  (28) 

yTk  =  [ru,«ik  cos(<t>'k  )  +  a!tk  cos(<|>k)]sin(>.k)  (29) 

zTk  =  fioci,  sin((J)'k  )  +  alt,  sin((j)k )  (30) 

F.  NONLINEAR  LEAST  SQUARES  ESTIMATION 

Since  the  (polar)  focal  plane  coordinates  of  the  observations  are  nonlinear  func- 
tions of  the  tactical  parameters,  a  nonlinear  least-squares  process  is  required.  This  is 
achieved  by  a  sequence  of  linear  least-squares  estimations.  Each  of  these  linear  least 
squares  estimations  is  based  on  a  one-term  (linear)  Taylor  series  expansion  of  the  obser- 
vations in  terms  of  the  tactical  parameters  A  general  requirement  of  ordinary  linear  least 
squares  is  that  the  noise  contaminating  the  observations  should  be  independent  and  have  a 
common  (preferably  Gaussian)  distribution 

In  this  case,  the  "polar"  focal  plane  observations  (3k,  tu)  are  transformed  into 
"Cartesian-like"  coordinates  (uk,  vk): 

uk  =-tan(rik)sin(Pk)  (31) 

vk  =-tan(r|k)cos(p\)  (32) 

so  that  the  coordinates  have  similar  behavior  with  respect  to  errors  and  noise  Either  r\  or 
sin(ri)  could  be  used  in  place  of  tan(r|).  but  since  tj<10°,  all  three  behave  similarly  The 
chosen  form  corresponds  to  a  gnomonic  projection  of  the  globe  onto  a  plane  tangent  at 
the  nadir  point    Figure  12  illustrates  this. 


21 


Figure  12.  Focal  Plane  Transformation 
Next,  it  is  necessary  to  determine  what  the  focal  plane  coordinates  would  be  for 
each  observation  if  the  TBM  were  actually  located  at  the  predicted  coordinates  given  in 
(28)  to  (30). 

The  expected  line-of-sight  is  computed: 

fAxA      (xT  -> 


Ay, 

^Azk. 


Urk  "zsk/ 


(33) 


In  order  to  determine  the  focal  plane  coordinates  of  azimuth,  P,  and  elevation,  r\, 
the  LOS  vector  given  in  (33)  must  be  rotated  into  the  observing  satellite's  local  vertical 
coordinate  frame,  (UEN): 


cos(ghak)      sin(ghak)     0YAxk^ 


-sin(ghak)    cos(ghak)    0 

V 


cos(5k)    0    sinCSJ"! 

0    1     0 
-sin(ok)    0    cos(Sk)A         0  0 

This  is  the  transpose  of  the  transformation  matrices  given  in  (7). 

Now  u  and  v  can  be  simply  expressed  in  (UEN)  terms  by  noting  that 

sin(Pk)  = 
cos(Pk )  - 


Aykj 

UzJ 


Ve7Tn7 

~Nk 

v^7+Nk2 


(34) 

(35) 
(36) 
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Expected  focal  plane  coordinates,  ( u  k  ,vk ) ,  are  then 

Ek 

u.  = 

Uk 

Nk 
v.  = 


(37) 


(38) 


Expected  focal  plane  coordinates  are  denoted  with  a  "  A  ".  Taking  the  difference  between 
actual  observations  and  expected  observations  generates  Su  and  8v  (found  on  the  left-hand 
sides  of  (41)  and  (42)): 

6uk  =  uk  -uk  (39) 

ovk=vk-vk  (40) 

Now  that  the  changes  in  the  focal  plane  coordinates  are  known,  they  can  be  related  to 

changes  in  the  tactical  parameters    The  one-term  Taylor  expansions  of  u  and  v  in  terms  of 

the  tactical  parameters  are: 

da 

da0 

dv 


6u 


du  „      da  __ 

5T0  - — 5Lh 

5T0      °     dL 


&J)()  + 8X,0  h —  5h0 

d$Q  dX0  dh0 


-5an 


(41) 


8v 


^     2T  ^V   21  ^     2X  (1W     21  tV     2U 

= 5T0  +  — 5L  + o<J)0  + o/i0  + oh0  + oa0 

cT0  cL  <3<t>0  f?X0  dh0  ca0 


(42) 


These  equations,  written  in  vector-matrix  format,  are  the  linear  least-squares 
model  that  forms  the  basis  of  the  tactical  parameter  estimation  process: 

fS"0 


fouk 

Uvk 
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duk 

duk 

^ik 

c\ik 

*0 

8L 

5Tfl 

a. 

a*« 

^0 

^h0 

aa0 

s*. 

^k 

avk 

avk 

dvk 

c-Vk 

^vk 

o^o 

l^To 

aL 

a*« 

ax0 

dh0 

da<J 

8h0 

(43) 


There  are  two  rows  of  the  center  matrix  in  (43)  for  each  of  the  k  =  (1  to  n)  observations 
for  a  total  of  2n  rows.  This  matrix  of  partials  is  denoted  the  "A"  matrix,  and  represents 
changes  in  the  focal  plane  coordinates  with  respect  to  changes  in  tactical  parameters 
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The  differences  in  the  focal  plane  coordinates,  5u  and  8v,  are  used  to  generate 
adjustments  to  the  tactical  parameters  by  solving  (43)  for  the  tactical  parameter  variations: 

6L 


84>,, 
ft, 

5h0 

6a0J 


^r'O 


(44) 


The  middle  term  in  (44)  is  the  pseudoinverse  of  A,  and  is  used  because  A  is  an  (2n  x  6) 
matrix,  and  is  generally  not  square  The  changes  to  the  tactical  parameters,  the  left-hand 
matrix  in  (44),  are  added  to  the  initial  estimate,  and  the  process  is  repeated  until  all  the 
components  of  the  left-hand  vector  in  (44)  are  sufficiently  small,  or  the  convergence  crite- 
rion is  met    The  result  is  the  correct  tactical  parameters  for  the  observed  TBM  launch 

The  matrix  of  partial  derivatives,  A,  is  the  key  to  the  whole  iteration  process  The 
development  of  this  matrix  involves  multiple  use  of  the  Chain  Rule  for  derivatives  The 
partial s  in  (4 1 )  can  be  expanded  to  yield: 

du  _  du  dty      du  dX      du  dh 
~'^Jl0^^x'^  +  ~dh~dT~0 


DT0 

du  _  du  dfy     du  dX     du  dh 

~dL~  ^~dL+~dX~dL+~di\^L 

du  __  du  d§      du  dk      du  dh 

aj>0    atj>aj>0    dxd$0    ahaj>0 

du   _du  6$      du  dX      du  dh 
^~~didk~0+~dk~dk^+a\~d)^ 

du  _  du  cfy      du  dX      du  dh 
dh0      cty  dh0      dX  dh0      oh  dh0 

du      du  d^      du  dX      du  <3h 

da0      cty  da0      dX  da0      cfa  da0 


(45) 
(46) 
(47) 
(48) 
(49) 
(50) 
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In  a  similar  manner,  the  partials  of  (42)  are  expanded: 

dv  _  dv  d§      dv  dk      dv  dh 

dv  _  dv  cfy      dv  dk     dv  dh 

dv  _  dv  ety      dv  dk      dv  dh 

dv   _  dv   dfy      dv  dk      dv  dh 

dv       dv  d^      dv  dk      dv  dh 
dh0  ~  a()  dh0      dk  dh0      dh  r?h0 

dv   _  dv  d$      dv  dk      dv   dh 
da0      dxj)  da0      dk  da0      dh  da0 

Patterns  and  structure  in  these  twelve  equations  are  more  easily  recognized 
these  equations  are  written  in  matrix  form: 


(51) 

(52) 

(53) 

(54) 

(55) 

(56) 
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(57) 
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db_ 

cTT0 

a* 

dL 
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dk 

dL 

dk 
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ax0 
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ah0 
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ah 
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a0 
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av 


av 
vah^ 


(58) 


vaa0y 

Note  here  that  (57)  and  (58)  have  identical  structure  The  left-hand  sides  are  the  deriva- 
tives of  the  focal  plane  coordinates  with  respect  to  the  tactical  parameters,  the  vectors  on 
the  right  are  the  derivatives  of  the  corresponding  focal  plane  coordinates  with  respect  to 
TBM  position  (latitude,  longitude,  and  height),  and  the  center  matrix  describes  the  deriva- 
tives of  the  position  with  respect  the  the  tactical  parameters  This  effectively  separates  the 
change  in  the  focal  plane  coordinates  with  respect  to  tactical  parameters  into  two  parts:  a 
matrix  that  describes  changes  in  position  due  to  trajectory,  and  a  vector  that  describes  the 
relation  between  focal  plane  observations  and  TBM  position 
The  matrix  is  given  a  name,  Ai: 
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aT0 
ah 


Ai  = 


a* 

dk 

ar0 

ar0 

a<|> 

dL 

dk 

dL 
dk 

atj> 

Wo 
dk 

ax0 

dk0 

a* 
ah0 

dk 

dh0 

d$ 

dk 

daQ 

daQ 

dL 

ah 

afro 
ah 

dk0 
ah 

ah 

dan 


(59) 
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so  thai  (57)  and  (58)  can  be  rewritten 


(  du> 

^Tn 

&l 

5L 

fda\ 

du 

a* 

du 

=  A, 

ax 

cfk. 

au 

du 

laJ 

dh0 

du 

V<3a0> 

(60) 


f  £v> 

c^n 

(5v 

a. 

ra/l 

_^ 

aj> 

^o 

-  A, 

rv 

0\ 

<"/. 

a0 

r.v 

GV 

^dh> 

ah0 

dv 

Vc5a0> 

(61) 


The  elements  of  A]  can  be  found  qualitatively     This  is  done  in  Appendix  A.  and 
the  results  are: 

7    cd  cos(a0)        dd     sin(a0) 


A, 


^0  feff 

1.5dP  cos(a0) 

refl 
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dsin(a0) 


cT0  reff  cos(<|)0) 
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r^COS^) 
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dcos(a0) 

reffcos((|>0) 


ah 

hp 
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(62) 
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The  values  of  the  elements  in  (62)  are  approximations,  but  approximate  magni- 
tudes and  correct  signs  are  all  that  are  required  for  the  iteration  process  to  converge  The 
benefit  of  making  this  simplification  for  Ai  is  the  reduction  in  computation  time  The 
complete  and  exact  formulae  can  be  extrapolated  from  Appendix  A 

Similarly,  vectors  on  the  right-hand  sides  of  (57)  and  (58)  can  be  separated  into 
simpler  components: 

du  _  du  dx     du  dy     du  dz 

^~^^  +  ^apaz"a* 

du  _  c\i  dx  dudy  du  dz 
dk      dx  dl.     dy  dk      dz  dk 

du  _  du  dx  du  dy  du  dz 
dh      dx  dh     dy  dh     dz  dh 

dv  _  dv  dx  dv  dy  dv  dz 
d§     dx  cty     dy  dfy     dz  dfy 

dv  dv  dx  dv  dy  dv  dz 
—  =  —  —  +  — - :r — ~~ 
dk      dx  dk     dy  dk     dz  dk 

dv  _  dv  dx     dv  dy     dv  dz 

dh      dx  dh     dy  dh     dz  dh 

As  before,  these  equations  can  be  written  in  matrix  form 
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(63) 
(64) 
(65) 
(66) 
(67) 
(68) 


(69) 


(70) 
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The  center  matrix  is  renamed  Av 


ax 

dy 

dz~ 

a* 
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dx 

dy 

dz 

ex 
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dx 

dy 

dz 

dh 

dh 

dh 

and  its  elements  can  be  approximated  by: 

-(reff  +  h)sin(<t>)cos(\)     -(reff  +  h)sin(<t>)sin(X)    (r^  +  h)cos(<t>) 
-(r^ +h)cos<<t>)sin(?i)     (r^  +  h)cos(4>)cos(X)  0 

cos(<j>)cos(A.)  cos(<t))sin(>t)  sin(<t>) 

The  vectors  on  the  right-hand  side  of  (69)  and  (70)  can  be  further  broken  down 
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The  center  matrices  in  (73)  and  (74)  are  identical,  and  renamed  A3 
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A3,  th  transformation  from  (XYZ)  to  (UEN)  has  been  derived  previously  in  (7): 


cos(ghak)    -sin(ghak)    0 

sin(ghak)      cos(ghak)     0 

0  0  1 


cos(6k)    0    -sin(8k) 

0  1  0 

sin(8k)     0     cos(6k) 


(71) 
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(73) 


(74) 
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The  final  step  is  to  determine  the  elements  in  the  vectors  on  the  right-hand 
(73)  and  (74)    These  are  easily  computed  from  (37)  and  (38): 

5u  __EL 

au~u2 

du_  __J_ 
dE  ~     U 

aN 

dv_  =  JN_ 

au  "  u2 

*  =  o 
dE 


£N  ' 
Now  (57)  and  (58)  can  be  rewritten: 
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1  5To 
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(84) 


voa0; 

Equation  (83)  is  used  to  compute  the  odd  rows  (k  =  odd)  of  A  in  (43),  and  (84)  is 
used  to  compute  the  even  rows  (k  =  even)  of  A  in  (43) 

Focal  plane  coordinates,  ((3k,  r|k).  are  transformed  into  Cartesian-like  focal  plane 
coordinates,  (uk,  vk)  (D.  (Circled  numbers  refer  to  steps  illustrated  in  Figure  13.)  Satellite 
positions  and  initial  estimates  of  tactical  parameters  are  used  to  generate  expected 
(theoretical)  focal  plane  coordinates,  (uk ,  vk )  3?Xa>,  and  the  A  matrix  X,  (which  is  denoted 
d(u,v) 


by 


in  Figure  13)  and  its  pseudoinverse  5      The  differences  of  the 


d(T0,L,<K,>.0,h0,a0) 

focal  plane  coordinates,  (5uk.  5vk),  (observed  minus  theoretical)  6  are  used  to  generate 
adjustments  to  the  tactical  parameter  estimates  by  computing  " 


(85) 


If  all  the  elements  of  the  left-hand  vector  (adjustments  to  the  tactical  parameters)  are  not 
sufficiently  small,  the  current  values  (j-1)  of  the  tactical  parameters  are  updated  ®  to  the 
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(86) 


If  the  iteration  process  has  converged,  the  iteration  process  is  terminated  and  the  current 
values  of  the  tactical  parameters  are  the  best  estimates  that  this  process  can  generate    This 
entire  process  may  be  visualized  by  a  flowchart  diagram  drawn  in  Figure  13. 
(M )  (gh*.  8,  R)         (To,  L.  4>0,  >.(!.  ho,  ao)^ 


(8 'I^  6L.  64>...  fi>.,..  6k,..  &0  - 


Figure  13    Flowchart  Diagram  of  Iteration  Process 

G.  BURNOUT  TIME  ESTIMATION 

The  estimation  of  burnout  time  is  based  on  the  TBM  profile  maximum  burn  time 
(W),  the  last  observation  time  (Tiaa),  and  the  next  potential  observation  time  (Tnex1),  had  it 
occurred    The  maximum  burn  time  according  to  the  profile  is: 

T_  =  To+t,^  (87) 

Two  cases  can  occur 
(a)ifTmax>Tne^then 

~  1 


\\ ; 


rO^-TL,) 


(88) 
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(b)ifTma,<Tnexlthen 


rfl^-T^) 


(89) 


A  figure  and  example  are  given  below  as  a  demonstration 
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Figure  14    Burnout  Time  Estimation 
For  the  sample  data  given  in  Table  1,  using  the  profile  given  in  Table  2.  t^  =  62  5 
seconds,  T)as  =  155  44,  and  the  next  potential  observation  is  TneTl  =  159  36     With  T0  = 
109  36,  the  maximum  burn  time  is  T^  =  171.86,  so  case  (a)  applies  and  burnout  time  is 
estimated  at  The  =  157.40. 

H.  STATE  VECTOR  GENTRATION 

The  state  vector  completely  defines  the  TBM's  position  and  velocity,  and  has  six 
elements  With  the  tactical  parameters  and  burnout  time  estimates,  generating  the  state 
vector  at  burnout  is  done  by  evaluating  the  position  and  velocity  equations  at  the  burnout 


=  T    -T 


(90) 
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Substitute  t*,  for  t  in  (1)  and  (2),  and  use  (3),  (4),  and  (22)  through  (25)  to  get  4h»,  foo, 

and  akbc,: 

dp0o  =ao  +a,tbo  +a2tj  +&,\J  +&4lJ 
h^=b.+b.t*+b2tJ+b3tfco3+b4tbo4 
^=(1-1.51^ 
hbo=d  +  L)hPte 


rcff 


<t>b0   =    --cos-'[cos(Gbo)sin((t>0)  +  sin(ebo)cos((|>0)cos(a0)] 


(1) 
(2) 
(3) 
(4) 

(22) 


(23) 


^,„  =  >.„  ^  sin- 


sin^  )sin(a0) 


(24) 
(25) 


cos(<|)bo) 
altbo  =  h0-hbo 

The  burnout  velocity  is  expressed  in  terms  of  speed  (Vbo),  flight  path  angle  (ybo), 
and  heading  (ctbo) 

vb(,  =  Vd2  +  h2  (91) 


7bo  =«an" 


»   h 


.     J  cos(<t>0)siiKao) 

abo  -  S,n  3  ^ 

v      cos(obo)      ; 
where  d  and  b  are  the  time  derivatives  of  d  and  h 
6  =  ™ 


a 


(92) 
(93) 

(94) 
(95) 
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(96) 


The  state  vector  at  burnout,  Tbo,  is: 

♦bo 

K 

K 

Ybo 

L  CALCULATION  OF  IMPACT  POSITION  AND  TIME 

{ The  equations  and  methods  in  this  section  are  adapted  directly  from  Chapter  6 
of  Bate,  et  al.)  [Ref  6]  The  ballistic  trajectory  is  modeled  in  three  phases  powered 
flight,  which  is  the  portion  that  DSP  observes,  free-flight,  which  is  a  portion  of  an  elliptical 
orbit,  and  re-entry,  which  is  the  portion  from  where  atmospheric  drag  becomes  significant 
until  missile  impact  The  powered-flight  phase  is  modeled  with  the  TBM  profile  polyno- 
mials presented  earlier  The  ellipse  traced  during  free-flight  is  simulated  using  inertial 
two-body  mechanics  Atmospheric  drag  effects  during  the  re-entry  phase  are  not  specifi- 
cally calculated  in  this  model,  but  are  somewhat  accounted  for  by  assuming  that  the  dis- 
tance traveled  over  the  earth  from  re-entry  to  impact  is  the  same  as  from  launch  to 
burnout  This  is  the  same  as  assuming  that  the  earth-central  angles  are  equal 

ere=eb0  (97) 

The  speed  of  the  TBM  during  re-entry  is  assumed  to  be  unaffected,  however  These 
approximations  are  recognized  artificialities,  but  are  simpler  than  calculating  ballistic 
coefficients  of  various  TBM's  and  modeling  atmospheric  density,  and  more  accurate  than 
assuming  the  missile  remains  on  its  elliptical  path  until  impact  Figure  15  illustrates  the 
geometry  and  some  of  the  quantities  used  in  the  equations 
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Figure  15    Ballistic  Trajectory  [Ref  6] 
The  equations  contained  in  Chapter  6  of  Ref  6  concern  ballistic  trajectories  and 
are  based  upon  inertia!  quantities  and  a  spherical  earth  model,  so  the  geocentric  position 
vector  and  inertial  speed  at  burnout  are  required     Geocentric  latitude  is  generated  from 
(26): 

fko  =  tan-,[(l-f)2tan(«bo)]  (26) 

The  earth's  radius  at  burnout  is  calculated  with  (27): 

r 'JlzSl an 

VO-f)2cos2(fbo)  +  sin2(<fr,bo) 
These  values  are  then  transformed  into  earth-centered,  Cartesian  coordinates  of  the  TBM 
xTto  =[rb«!bocos((t),bo)  +  altbocos(<|)bo)]cos(?.bo)  (28) 

y^  =[rtoe.lboCOs(fbo)  +  altbocos(«J)bo)]sin(Xbo)  (29) 

ztw  =  W  sin((^'b0)  +  altb0  sin((t>bJ  (30) 
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which  define  the  geocentric  TBM  position  at  burnout: 

In  an  inertial  reference  frame,  the  TBM  has  an  initial  eastward  velocity  at  launch 
due  to  the  earth's  rotation,  v0    This  velocity  can  be  expressed  as: 

v„  =  rloc.lboO>ecos((J>;o)  (98) 

where  (Oe  =  15°/hour,  the  rotation  rate  of  the  earth  The  burnout  velocity  can  be  broken 
down  into  three  inertial  (UEN)  components  In  the  inertial  frame,  the  initial  eastward 
velocity  is  added  to  the  eastward  component: 

Vu=Vbosin(ybo)  (99) 

VE  =  Vbocos(Y1Jsin(alJ  +  v0  (100) 

VN  =  Vbocos(Ybo)cos(ab<))  (101) 

Inertial  quantities  are  denoted  with  a  bar.  "  _  ".  The  magnitude  of  this  vector  is  the  iner- 
tial speed,  Vb(, 

V^VV-'vV+'vy  (102) 

The  inertial  flight  path  angle  and  heading  may  now  be  calculated: 


(103) 


(  V  ^ 
a^tan'^  (104) 

Now  that  the  inertial  quantities  are  known,  a  non-dimensional  parameter,  Q,  is  defined  as 
the  squared  ratio  of  the  inertial  speed  of  the  TBM  to  circular  orbit  speed  at  that  position: 

V  2r 
Q=     *    te  (105) 

P- 

where  u  =  398601  2  — -   the  gravitational  parameter  for  the  earth,  and  rbo  =  irb*|,  the 
sec 

magnitude  of  the  TBM  position  vector 
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The  eccentricity  of  the  ballistic  "orbit"  is: 

e^l+QCQ^cos2^)  (106) 

The  free-flight  "earth  central  angle",  *F,  is  defined  by  Bate,  et  al ,  as  that  earth- 
central  angle  that  the  missile  traverses  between  burnout  and  re-entry  (see  Figure  1 5)  This 
definition  is  modified  by  adding  the  angle  traversed  during  reentry,  0re  to  the  BMW  de- 
fined angle,  *¥  This  implements  the  assumption  that  the  TBM's  trajectory  is  modified  by 
drag    The  new  definition  for  the  angle,  4/,  is: 

y  =  a,ot"(lzQ£ggi8)W  (107) 

The  eccentric  anomaly,  E,  and  the  semi-major  axis,  a,  of  the  ballistic  trajectory  are: 


E  =  cos  ' 


,   e  -  cos 

V2 


1  -  ecos  — 


(108) 


a  =  — *-  (109) 

2-Q 

Assuming  the  speed  of  the  TBM  is  unafTected  by  drag  during  re-entry,  the  time  of 

free  flight  (burnout  to  impact)  is  defined  as 

tff=2  ;— [rc-E  +  esin(E)]  (110) 

the  time  of  impact  from  launch  is 

tu„  =  tbo  +  tff  (111) 

and  the  time  of  day  of  impact  is: 

Tta«T0  +  U  (112) 

Using  the  Law  of  Cosines  from  spherical  trigonometry,  latitude  at  impact,  $'„„,  is 

4>L  =  sin,[sin(<t)^)cos(^)- 005(271-51^)]  (113) 

This  can  be  transformed  back  into  the  WGS-84  latitude: 
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This  method  is  not  an  exact  transformation,  but  it  is  a  close  approximation  that  does  not 
affect  the  error  analysis  A  better  model  should  be  used  if  more  precise  latitude  of  impact 
is  desired 

Using  the  Law  of  Cosines  again,  and  after  some  algebra,  the  longitude  traversed, 
AX,  is: 

^^■.(^-■O^)].  (]]5) 

I        cos(Ocos(*L.)       J 
So,  the  impact  longitude,  ^lun,  is: 

Xun  =  Xbo+AX  (116) 

Better  models  for  the  ballistic  trajectory  could  possibly  be  employed     The  goal  of 

this  analysis  is  not  to  precisely  determine  the  actual  impact  zone,  but  to  determine  the 

effects  of  error  sources  upon  the  result    Thus,  the  effects  found  using  this  model  should 

be  applicable  to  other  ballistic  trajectory  models 
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IV.  ERROR  SOURCES 

Noise  (error)  is  present  in  every  aspect  of  the  algorithm  presented  in  Chapter  3 
Every  measurement,  constant,  coefficient,  model,  equation,  approximation,  and 
assumption  impart  some  finite  quantity  of  uncertainty  upon  the  result  Even  if  a  quantity  is 
exact  to  its  limit  of  accuracy,  the  simple  fact  that  it  is  represented  by  a  finite  number  of 
digits  limits  the  precision  of  that  quantity  For  example,  if  a  ruler  has  one-sixteenth  inch 
gradations,  the  accuracy  of  any  measurements  made  with  the  ruler  is  ±one-thirty-second  of 
an  inch  This  means  that  there  are  no  "perfect"  values  in  the  algorithm,  and  each  value 
used  is  a  source  of  uncertainty,  termed  an  error  source 

Each  error  source  usually  has  one  or  more  underlying  causes  A  complete  error 
analysis  would  break  each  error  source  down  to  its  fundamental  level,  and  model  each 
level  correctly  An  analogy  for  this  is  the  layers  in  an  onion  The  onion  can  be  viewed 
externally  as  a  whole,  but  can  also  be  peeled,  layer  by  layer,  to  reveal  deeper,  underlying 
matter  that  support  the  external  skin  In  this  section,  some  of  the  underlying  causes  of  the 
overall  error  sources  are  identified,  but  in  these  analyses,  only  the  overall  error  magnitudes 
will  be  modeled  and  analyzed 

The  obvious  origin  of  any  errors  present  is  the  observational  data  since  it  is  the 
starting  point  for  the  algorithm  These  data  are  physical  measurements  of  three  types 
time,  focal  plane  measurements,  and  satellite  position  and  orientation. 

Time  errors  can  be  caused  simply  by  having  more  than  one  clock  being  referenced 
as  a  source  of  time  measurement  Imperfect  synchronization  between  clocks'  time  and 
time  passage  rate  are  obvious  error  sources  Time  delays  caused  by  radio  transmission  of 
data  due  to  distance,  atmospheric  refraction,  and  relative  motion  Doppler  effects  may  add 
another  time  error  In  the  extreme,  the  source  of  time  error  can  be  traced  all  the  way  back 
to  our  measurement  of  time  passage  with  respect  to  the  celestial  sphere,  which  is  not 
fixed  The  magnitude  of  time  errors  is  relatively  small,  and  getting  smaller  as  time 
measurement  improvements  are  being  made  continually 


The  largest  source  of  time  error  enters  the  algorithm  at  the  burnout  time  estimation 
phase  The  10-second  sampling  rate  limits  the  accuracy  of  the  estimate  to  ±  5  seconds  for 
single  satellite  observations  The  accuracy  of  the  profile's  value  for  t^  also  comes  into 
play  in  the  estimate  The  burnout  time  is  crucial  to  determining  the  ballistic  trajectory, 
since  the  TBM  is  undergoing  its  greatest  acceleration  changes  during  the  final  seconds  of 
powered  flight  Small  errors  in  the  estimate  of  burnout  time  thus  result  in  larger  errors  in 
the  estimate  of  impact  time  and  position 

LOS  measurement  errors  can  arise  from  many  noise  sources,  mainly  of  two  types: 
attitude  uncertainties  and  IR  radiation  measurement  errors  The  attitude  of  the  spacecraft 
needs  to  be  precisely  known  (and  it  simplifies  the  process  if  it  is  very  stable)  A  28 
microradian  error  in  pitch  or  roll  pointing  accuracy  of  a  geosynchronous  satellite  equates 
to  a  1  kilometer  error  on  the  surface  of  the  earth,  at  least  35786  kilometers  away  Any 
attitude  control  system  inaccuracy  effects  are  amplified  by  the  geosynchronous  altitude 
Vibrational  noise  from  spinning  momentum  wheels,  uncertain  knowledge  of  the  mass 
properties  of  the  yaw-spinning  satellite,  star  catalog  and  star  sensor  measurement  errors, 
thruster  misalignments  and  disturbance  torques,  tachometer  errors,  and  mechanical 
misalignments  all  contribute  to  the  total  attitude  uncertainty  The  knowledge  of  the 
telescope  boresight  axis  alignment  with  the  satellite's  reference  frame  is  defined  by  the 
design  and  manufacturing  of  the  DSP  satellite,  and  changes  slightly  with  thermal 
variations 

The  IR  radiation  measurements  of  intensity  and  angle  of  arrival  (AOA)  also  have 
several  underlying  error  sources  Locations  of  the  individual  photoelectric  cells  on  the 
focal  plane  array  are  recorded  in  what  is  termed  the  "Focal  Plane  Vector  Table"  (FPVT) 
The  positions  of  the  PEC's  change  as  the  satellite  heats  and  cools  with  varying  sun- 
satellite  orientations,  causing  a  warping  of  the  focal  plane,  but  the  FPVT  does  not  account 
for  the  changes  real-time  In  addition,  detector  noise,  refraction,  and  IR  attenuation  due 
to  clouds  and  water  vapor  all  add  to  the  total  measurement  error. 

Satellite  position  measurements  from  AFSCN  are  used  to  update  orbit  ephemeris 
data  The  measurements  are,  of  course,  inexact,  but  those  errors  are  compounded  over 
time  because  the  updates  occur  only  once  weekly    The  ephemeris  data  is  propagated  to 
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determine  satellite  position  during  the  week  between  updates  DSP  is  a  geosynchronous 
satellite,  so  hs  orbit  does  not  change  greatly  in  that  period,  but  errors  in  the  propagation 
model  can  cause  the  predicted  position  to  be  different  from  reality 

Keeping  all  of  these  underlying  components  in  mind,  the  errors  are  modeled  in  the 
tactical  parameter  estimation  algorithm  The  real-world  error  component  magnitudes,  bias 
and  random  distributions  may  be  different  from  those  simulated,  but  the  overall  effects 
manifest  themselves  in  a  manner  close  to  the  error  models  It  should  be  possible  to 
extrapolate  the  results  obtained  in  this  study  to  different  errors  by  properly  scaling  the 
different  magnitudes  and  distributions  It  must  be  emphasized  that  since  the  goal  of  this 
study  is  not  to  exactly  determine  the  tactical  parameters  at  launch,  state  vector  at  burnout, 
or  impact  time  and  position  The  purpose  is  to  determine  the  contribution  of  the  error 
sources  upon  these  values,  and  this  can  be  done  with  approximate  models  for  the  errors 
If  the  intended  goal  is  to  calibrate  the  system  to  eliminate  bias  errors,  then  the  sources  of 
the  errors  must  be  determined  more  exactly  to  determine  what  is  observable  (measurable) 
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V.  NUMERICAL  ANALYSIS 

The  algorithm  described  in  Chapter  III  was  programmed  into  MATLABIM  Three 
error  sources  (time,  satellite  position,  and  LOS)  were  simulated  and  added  to  the 
observational  data,  first  separately,  and  then  combined.  The  effects  of  the  errors  upon  the 
results  are  analyzed  at  three  points:  launch  position,  burnout  position,  and  impact 
position  Five  Middle  Eastern  capital  cities  were  chosen  as  fictitious  launch  sites  to 
determine  the  effects  of  TBM  launch  position  upon  the  accuracy  of  the  results. 

The  model  for  time  error  is  a  random  uniform  distribution  between  0  and  1 
milliseconds.  For  satellite  position,  a  random  normal  distribution  with  zero  mean  and 
standard  deviation  of  200  meters  was  added  to  each  of  the  components  (R,  gha,  5).  This 
is  approximately  a  one  standard  deviation  sphere  with  radius  346  meters  The  LOS  error 
was  modeled  as  a  random  normal  distribution  with  a  standard  deviation  of  5  microradians, 
added  to  both  (3  and  r\  components  Histograms  of  sample  error  distributions  are  shown 
in  the  following  three  figures 


Figure  16.  Histogram  of  Time  Error  -  Uniform  Distribution 
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Figure  17    Histogram  of  Satellite  Position  Error  -  Normal  Distribution 


-1         -0.5  0  0.5  1 

Focal  Plane  Error  (radians) 


Figure  18    Histogram  of  LOS  Error  -  Normal  Distribution 


46 


Errors  added  to  the  observations  had  the  effect  of  creating  an  ellipse  of  incorrect 
data  points  around  the  true  position,  termed  an  "error  ellipse"  The  1  -<r  ellipses  shown  in 
the  next  figures  are  determined  using  statistical  methods  [Ref  7]  For  the  launch  ellipse,  a 
matrix  of  two  columns  was  collected  during  the  simulation: 

V.      4*0.. 


lpts  = 


O-ro-o 


(117) 


where  lpts  =  launch  points,  the  matrix  of  launch  longitudes,  Ao,  and  latitude,  <t>0.  The  first 
row  is  the  generated  true  launch  position,  the  second  row  is  the  calculated  launch  position 
with  no  error  added  to  the  observational  data,  and  the  remaining  999  rows  are  the 
calculated  launch  positions  with  errors  added  to  the  observational  data  The  first  two 
rows  are  error-free,  and  are  removed  to  produce  a  matrix  of  only  "corrupt"  positions  To 
determine  the  ellipse  dimensions,  a  and  b  (semi-major  and  semi-minor  axes  lengths),  the 
eigenvectors  and  values  of  the  covariance  of  this  matrix  were  calculated  Twice  the 
square  root  of  the  eigenvalue  diagonal  elements  produces  a  and  b: 


=  2' 


v/eigenvaiue(l,l)  0 

0  <Jeigenva\ue(2,2) 


(118) 


and  the  eigenvector  is  the  direction  cosine  matrix  for  rotating  the  ellipses  from  the  primary 
axes  The  impact  ellipses  are  calculated  similarly,  with  the  initial  data  matrix  being 
composed  of  impact  longitudes  and  latitudes,  instead  of  launch  The  ellipses  can  then  be 
plotted  using  the  ellipse  equation: 


y_ 

b2 


(119) 


and  letting  x  vary  between  ±a,  and  letting  y  be  the  dependent  variable.  Multiplying  the 
resultant  (x,  y)  coordinates  by  the  eigenvector  matrix  rotates  them  to  the  correct 
orientations   The  semi-axes  are  in  units  of  radians,  since  the  positional  values  are  in 
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radians.   Thus,  the  (x,  y)  plots  produced  are  the  longitude  and  latitude  of  the  points,  once 
the  values  are  converted  to  degrees 

For  the  burnout  position,  the  same  method  is  used,  but  the  data  point  matrix  has 
three  columns,  the  third  being  the  altitude  of  the  TBM  at  burnout.  The  covariance  method 
works  as  well  in  three  dimensions,  producing  a,  b,  and  c,  the  three  semi-axes  for  the 
ellipsoid  formed  The  coordinates  are  calculated  in  the  same  manner,  using  the  equation 
for  an  ellipsoid 


y_ 

b2 


:1 


(120) 


and,  again,  rotating  by  pre-multiplying  by  the  eigenvector  matrix 

The  following  figures  are  representative  of  all  the  simulated  cases,  but  this 
particular  case  is  a  300-km  TBM  launched  from  Baghdad  heading  0°,  with  the  combined 
erorrors 

Launch  Position 
33.337  ■ 


44.43 


44.432         44434         44  436 
Longitude  (degrees) 


Figure  19    Launch  Position  Error  Ellipse 
The  ellipse  drawn  on  top  of  the  data  points  encompasses  the  data  points  within  one 
standard  deviation  of  the  mean,  assuming  the  distribution  is  normal.    Since  the  time  error 
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has  a  uniform  distribution,  l-o  is  not  the  expected  66%  for  time  error  cases  This  method 
is  used  anyway  for  lack  of  a  better  one  and  to  maintain  a  common  baseline  for 
comparison.  As  it  turns  out,  the  time  error  contribution  is  very  small,  and  thus  does  not 
effect  the  combined  error  case  distribution. 

The  position  at  burnout  is  plotted  in  three  dimensions,  using  altitude  as  the  third 
coordinate  to  form  an  ellipsoidal  error  volume 

Position  at  Burnout 


Latitude 


Longitude  (degrees) 


Figure  20    Position  at  Burnout 
The  three  ellipses  show  the  outline  of  the  ellipsoidal  volume     It  is  difficult  to 
perceive  the  ellipse  orientations,  but  they  are  mutually  orthogonal,  (necessarily,  since  they 
correspond  to  the  eigenvectors  of  distinct  eigenvalues) 
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The  impact  ellipse  is  similar  but  larger  than  the  launch  position  ellipse,   due  to  the 
propagation  of  velocity  errors  from  burnout  to  impact  and  the  error  in  the  burnout  time 

Impact  Position 
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Figure  2 1    Impact  Position  Error  Ellipse 

The  three  previous  figures  illustrate  the  error  ellipses(oids),  but  showing  all  the 
data  obtained  in  this  format  is  unwieldy.  The  sizes  (areas  and  volumes)  of  the  positional 
errors  are  the  important  quantities  for  identifying  which  error  source  has  the  largest  effect 
Thus,  the  data  have  been  processed  to  determine  these  quantities,  and  are  plotted  next. 

The  data  for  the  300-km  TBM  was  collated  by  city,  heading,  and  error  source. 
The  ellipse(oid)  dimensions  calculated  must  be  equated  to  the  actual  distance  on  the 
ground  to  determine  the  area  of  the  ellipses  To  do  this,  the  angles  must  be  in  radians,  and 
are  multiplied  by  riocai  The  distance  between  lines  of  longitude  shrinks  as  the  latitude 
moves  away  from  the  equator,  so  the  longitudinal  distance  must  be  multiplied  by  the 
cosine  of  the  geocentric  latitude,  (J)'    Thus  the  equation  for  ellipse  area  becomes: 

area  =  7irlocal2abcos((|)')  (121) 
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The  ellipsoid  volume  is  not  as  simple,  so  a  parallelepiped  with  dimensions,  a  x  b  x 
c,  is  used  to  approximate  it: 

volume  =  abc  (122) 

These  areas  and  volumes  are  plotted  in  the  following  figures,  separated  by  error 
source  and  heading    Each  launch  point  is  represented  by  a  different  color: 
Aden  =  blue 
Baghdad  =  black 
Damascus  =  green 
Riyadh  =  red 
Tehran  =  cyan 
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Figure  22    Time  Error  Effects  Upon  Launch  Ellipse  Area 
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Figure  23    Time  Error  Effects  Upon  Burnout  Ellipsoid  Volume 
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Figure  24    Time  Error  Effects  Upon  Impact  Ellipse  Area 


52 


0  50  100  150  200  250  300  350 

Heading  (degrees) 

Figure  25    Satellite  Position  Error  Effects  Upon  Launch  Ellipse  Area 
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Figure  26    Satellite  Position  Error  Effects  Upon  Burnout  Ellipsoid  Volume 
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Figure  27    Satellite  Position  Error  Effects  Upon  Impact  Ellipse  Area 
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Figure  28    LOS  Error  Effects  Upon  Launch  Ellipse  Area 
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Figure  29    LOS  Error  Effects  Upon  Burnout  Ellipsoid  Volume 
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Figure  30    LOS  Error  Effects  Upon  Impact  Ellipse  Area 
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Figure  3 1    Combined  Error  Effects  Upon  Launch  Ellipse  Area 
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Figure  32.  Combined  Error  Effects  Upon  Burnout  Ellipsoid  Volume 
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Figure  33  Combined  Error  Effects  Upon  Impact  Error  Ellipse 
It  is  immediately  apparent  from  these  plots  that  the  launch  site  has  little  effect  upon 
the  size  of  the  error  ellipses  and  ellipsoids  It  is  also  obvious  that  the  size  is  dependent 
upon  TBM  heading  and  observation  geometry  The  relative  magnitudes  of  the  error 
effects  are  not  as  obvious,  so  a  plot  of  four  error  ellipses  is  shown  in  the  next  figure,  with 
each  ellipse  formed  by  a  different  error  source 
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Figure  34.  Launch  Error  Ellipses  for  Each  Error  Source 
From  this  plot  it  is  easily  recognized  that  the  errors  can  be  ranked  in  order  of 
effect: 

1 .  focal  plane  errors 

2  satellite  position  errors 

3  time 

All  the  error  sources  together  produce  the  largest  ellipse,  as  would  be  expected  from  the 
superposition  principle 

This  is  further  illustrated  by  Table  3  Each  error  source  is  a  row:  "time"  is  the 
time  error  runs,  "satpos"  is  the  satellite  position  error  runs,  "LOS"  is  line  of  sight,  and 
"combin ."  is  the  combined  errors  runs 


ERROR 

LAUNCH  ELLIPSE  AREA 
(km2) 

BURNOUT  ELLIPSOID 
VOLUME  (km3) 

IMPACT  ELLIPSE  AREA 
(km2) 

nun 

mean 

max 

min 

mean 

max 

min 

mean 

max 

time 

3.26e-8 

2.59e-7 

1.22e-6 

2.9e-ll 

3.6e-10 

1.33e-9 

1.57e-6 

2.61e-5 

1.41e-4 

satpos 

3.77e-3 

5.27e-3 

8.68e-3 

6.04e-4 

1.14e-3 

2.32e-3 

6.12e-3 

7.92e-2 

3.27e-l 

LOS 

2.97e-2 

8.45e-2 

1.83c- 1 

2.37e-2 

9.33e-2 

2.09e-l 

3.05eO 

1.02el 

1.95el 

combin. 

3.91e-2 

1.01e-l 

2.97e-l 

4.44e-2 

1.25e-l 

4.50e-l 

3.35e() 

1.38el 

4.06e  1 

Table  3    Error  Ellipse  Areas  and  Ellipsoid  Volumes 
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This  is  also  shown  graphically  in  the  following  plots: 
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Figure  35    Launch  Ellipse  Area  Comparison 
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Figure  36.   Burnout  Ellipsoid  Volume  Comparison 
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Figure  37    Impact  Ellipse  Area  Comparison 
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The  sum  of  the  error  sources'  positional  error  sizes  is  less  than  the  combined  error 
source's  positional  error  size  This  implies  a  nonlinear  interaction  of  error  effects  with 
each  other  This  is  a  reasonable  result  considering  the  complexity  of  the  simulated  system, 
and  noting  that  most  real  systems  exhibit  nonlinear  behavior. 

Taking  advantage  of  the  ease  of  modifying  MATLAB™  code,  various  changes  to 
the  present  DSP  system  model  were  put  into  simulation  to  determine  the  effects  upon  the 
accuracy  of  the  results  For  every  case,  Baghdad  is  the  launch  site  and  the  observational 
data  were  modified  by  combined  error  sources 

The  first  case  is  the  control  case  with  nominal  system  parameters  It  is  listed  as 
"Bagcom"  to  represent  "Baghdad  combined  errors",  and  is  used  as  a  baseline  case  for 
comparison  The  second  case,  "Synchr "  shows  the  effects  of  synchronizing  the  spins  of 
the  satellites  so  that  the  satellites  scan  a  single  area  of  interest  within  one  second  of  each 
other  For  the  third  case,  "Molniya",  the  satellite  at  70°  GHA  was  elevated  to  63  4°  decli- 
nation, to  simulate  a  Molniya  +  geostationary  viewing  geometry  The  remaining  cases 
were  simulating  faster  scan  rates  from  10  seconds  down  to  2.5  seconds  The  numerical 
results  are  shown  in  Table  4. 


ERROR 

LAUNCH  ELLIPSE  AREA 
(km2) 

BURNOUT  ELLIPSOID 

VOLUME  (km3) 

IMPACT  ELLIPSE  AREA 
(km2) 

nun 

mean 

max 

mm 

max 

max 

nun 

mean 

max 

Bagcom 

4.2U-2 

l.Ole-1 

2.74c-] 

4.46e-2 

1  09e-l 

2.60c-] 

3.83eO 

1.40el 

4.04el 

Svnchr 

440e-2 

7.07e-2 

1.10e-l 

2.09e-2 

4.89e-2 

8.61e-2 

3.94e() 

7.82e0 

1.52e0 

Molruva 

4.39e-2 

7.58e-2 

1.25e-l 

261e-2 

5.57e-2 

104e-l 

2.91e0 

7.96e0 

1  85el 

lOssc 

4.21e-2 

1.01e-l 

2.74e-l 

4.46e-2 

1.09e-l 

2.60e-l 

3.83eO 

1.40el 

4  04el 

7.5  ssc 

2.62e-2 

7.44e-2 

1  95e-l 

2.53e-2 

4.55e-2 

1.09e-l 

403e0 

8.72e0 

2.08el 

5ssc 

2.48e-2 

5.02e-2 

1  52e-l 

1.23e-2 

2.52e-2 

6.52e-2 

252eO 

5.84c0 

1.64cl 

2.5ssc 

1.18e-2 

2.36e-2 

3.48e-2 

5.05e-3 

9.21e-3 

1  31e-2 

1.03e0 

2.95e0 

4.58e0 

Table  4    Various  Case  Comparison 
The  effects  of  the  modifications  are  more  easily  observed  graphically 
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Figure  38  Comparison  of  Mean  Launch  Ellipse  Areas 
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Figure  39    Comparison  of  Mean  Burnout  Ellipsoid  Volumes 
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Figure  40    Comparison  of  Mean  Impact  Ellipse  Areas 

Both  modifications  have  the  effect  of  decreasing  the  areas  and  volume  by  about  a 
third  This  is  the  expected  result  for  the  Molniya  case,  since  the  viewing  geometry  is  more 
three-dimensional  with  two  satellites  on  the  equator  and  the  middle  satellite  near  Molniyan 
apogee. 

The  synchronized  spin  results  were  a  surprise,  however.  The  effect  upon  the 
launch  ellipse  area  is  fine,  but  since  the  burnout  time  estimation  should  have  been  less 
accurate,  the  burnout  ellipsoid  volume  and  impact  ellipse  area  were  expected  to  remain 
unchanged  at  best.  Obtaining  simultaneous  observations  would  have  obvious  advantages, 
since  they  could  be  processed  differently  and  be  used  to  determine  the  position  exactly  for 
discrete  instants  in  time  during  the  boost  trajectory  However,  these  observations  were 
not  processed  differently,  and  were  not  simultaneous  ~  only  close  together  (within  one 
second  of  each  other).  Still,  the  result  may  be  a  manifestation  of  some  mathematical 
process  that  causes  the  increase  in  accuracy. 

The  following  plots  show  the  effects  of  increasing  the  scan  rate,  or  alternmely, 
decreasing  the  time  between  observations.  This  allows  more  observations  to  be  made 
during  the  boost  phase 
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Figure  42.  Scan  Rate  Effects  on  Burnout  Ellipsoid  Volume 
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Figure  43.   Scan  Rate  Effects  on  Impact  Ellipse  Area 
The  effects  seem  to  be  linear  decreases  upon  the  launch  and  impact  ellipse  areas, 
and  quadratic  decreases  on  the  burnout  ellipsoids    This  result  is  in  perfect  agreement  with 
intuitive  predictions. 

In  summary,  the  numerical  results  are  generally  in  accordance  with  expected  re- 
sults. As  a  qualitative  check  of  the  quantitative  results,  the  next  chapter  is  devoted  to  a 
qualitative  analysis  of  a  simplified  case  for  comparison. 
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VI.  QUALITATIVE  ANALYSIS 

Starting  with  the  equations  presented  in  Chapter  3,  a  qualitative  analysis  has  been 
done  to  determine  the  expected  trends  in  the  numerical  analysis  The  problem  is  simplified 
by  assuming  that  the  equations  are  exact  and  that  there  is  only  one  satellite  The  effects  of 
time  errors  are  not  analyzed,  so  the  problem  is  to  determine  how  satellite  position  and 
focal  plane  error  effects  should  effect  the  results 

TBM  geocentric  longitude  and  latitude  can  be  expressed  in  (XYZ)  coordinates 

X  =  tan   'li  (11) 


These  (XYZ)  coordinates  can  be  expressed  in  terms  of  satellite  position  and  focal  plane 
coordinates  using  the  LOS  projection  equation  (8) 

yT    =Rt+pe  (8) 

Vz,J 

Rewriting  the  terms  on  the  left-hand  side  in  terms  of  R,  gha,  5.  r\,  and  p  generates  a  rather 
lengthy  equation,  and  is  given  in  Appendix  B  The  next  task  is  to  find  the  partial 
derivatives  of  the  TBM  position  with  respect  to  the  five  error  terms,  using  the  Chain  Rule 
Taking  the  partial  with  respect  to  R.  for  example 

c5R~axc5R  +  ayaR  +  azaR 

dR      d\  8R  +  ay  dR      cz  dR 
Again,  expression  this  equation  written  explicitly  in  terms  of  the  error  sources  becomes 
tedious,  and  can  be  found  in  Appendix  B    The  results,  however,  are  interesting  and  are 
shown  as  plots    The  horizontal  axis  pointing  left  is  r\.  the  horizontal  axis  pointing  right  is 
p,  and  the  vertical  axis  is  the  differential 
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The  partial  derivatives  with  respect  to  gha  turn  out  to  be  simple 

ST-1  <125> 

dgha 

lr  =  0  (126) 

dgha 

The  observations  one  can  make  from  analyzing  these  plots  is  that  most  of  the 
partial  derivative  plots  are  flat  and  close  to  zero  When  ^approaches  8  6°,  however,  some 
of  the  values  begin  to  get  large  This  happens  when  the  observed  IR  event  is  close  to  the 
limb  of  the  earth  as  viewed  from  a  geosynchronous  satellite  There  also  appear  to  be 
transitions  when  r\  approaches  8  6°  and  b  is  0°  or  180°. 

Figures  44  and  45  show  that  satellite  radius  error  effects  are  of  such  small 
magnitude  that  they  are  neglible  Figures  46  and  47  show  that  satellite  declination  error 
effects  are  more  important  to  determine  latitude  than  longitude,  except  when  r\ 
approaches  8  6°  Figures  48  and  49  show  that  (3  LOS  errors  affect  latitude  determination 
more  than  longitude,  except,  again,  when  r|  approaches  8  6°  Figures  50  and  51  show  that 
r\  LOS  errors  have  the  largest  effect  of  all,  and  the  partials  exhibit  the  same  behavior  near 
the  limb 

These  plots  of  formulae  give  an  indication  of  what  magnitude  the  error  ellipse 
should  have,  given  the  magnitude  of  the  error  is  known,  by  using  that  magnitude  in 
equations  similar  to  (127) 

^^J****^**^  (127, 

dR  Vex  dR     dy  dR     dz  dRJ 

where  R  can  be  any  of  the  five  error  quantites  (R,  6,  gha,  r\,  or  p) 

To  test  this  assumption,  substitute  the  magnitude  of  the  errors  used  in  the 
simulation,  and  compare  the  result  with  the  numerical  analysis    The  errors  are: 

AR  =  200m  (128) 

.(    0.200  km    "\  „„rtx 

A6  =  tan"'  (129 

U2164.17kmJ 
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,     0.200  km    ^  „„„, 

Agha  =  tan' (130 

U2164.17kmJ 

Ar|  =  5  ujadians  (131) 

A(3  =  5  uxadians  (132) 

The  generalized  equations  for  the  partial  derivatives  of  the  geocentric  longitude 

and  latitude  can  be  used  to  find  the  derivatives  for  a  specific  location  by  using  specific 

focal  plane  coordinates    For  Baghdad,  viewed  from  the  satellite  at  70°  gha,  in  particular, 

these  coordinates  would  be: 

r)  =  01216  radians  (133) 

P  =  38546  radians  (134) 

Using  these  values,  the  partial  derivatives  with  respect  to  the  five  error  sources 

become 

(135) 
(136) 
(137) 
(138) 
(139) 
(140) 
(141) 
(142) 
(143) 


a 

^^-.^  rad 

-0.573 

ap~ 

rad 

at>' 

~m~m     rad 

=  -0.753 

ap = 

rad 

ex 

rad 

-13.788 

<?n 

rad 

a<t>'  _ 

r  r.,  rad 

=  5.96 

dn  = 

rad 

ex 

„    rad 

0.656 

58  ~ 

rad 

cx|>' 

~    .„  rad 

=  -0.658 

"as"" 

rad 

dX 

i  rad 

-3.663  xl  0_l 

~dR~ 

m 

a4>' 

rad 

=  1.583  xl  0"' 

ai " 

m 

a 

rad 

cgha 

rad 

TT  =  0^  <144> 

ogha        rad 
Multiplying  (135)  through  (144)  by  the  appropriate  error  quantities  produces  the 
geocentric  longitude  and  latitude  error  values 

A<t^  =  3.166x10^  rad  (145) 

A^.  =0rad  (146) 

A<fr;  =-3.121x10-*  rad  (147) 

A<J>;  =2.98x10-'  rad  (148) 

A<t>p  =  -3.765x10^  rad  (149) 

AXR  =  -7.326x10^  rad  (150) 

A^gh.  =  4.743x10^  rad  (151) 

AXS  =  -3.112  xlO"*  rad  (152) 

AX „  =-6.894x10  ?  rad  (153) 

AXP  =  -2.865x10^  rad  (154) 

Grouping  these  in  terms  of  satellite  position  and  focal  plane  and  summing  the  absolute 
values 

^Luuucpo^ooe™  =  W*  +  W*.  +  A(f>;  =  6.287  x  \0~*  rad  (155) 

&**»,**.-*  =  AXR  +  AXgh.  +  AX6  =  1518x10  ?  rad  (156) 

A*focdpu..«or  =  A<K  +  a4>p  =  3-357  x  10  ?  rad  (157) 

A>ifoe.lpUaeCTIOT  =  AX,  +  AXp  =  7.181  x  10 -'  rad  (158) 

These  values  are  in  radians,  which  can  be  equated  to  kilometers  on  the  surface  of 

the  earth  by  multiplying  latitude  by  the  earth's  radius  and  longitude  by  the  earth's  radius 

and  the  cosine  of  the  latitude    For  this  analysis,  earth's  radius  is  assumed  to  be  6378.137 

kilometers 

A4>Lntepo.ao.««  =  (6287  x  10^  rad)  6378.137  km  =  4.010  x  10"2  km  (159) 

A^uncpcmoccTo,  =  (1518  x  10"3  rad)cos(33.333°)6378.137  km  =  8.089  x  10  2  km  (160) 
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A^ocpta..™.  =  (3357  x  10  ?  rad)  6378  137  km  =  2  141  x  10  '  km  (161) 

AXf0C-pU.eem„  =  (7.181  x  10 -'  rad)  cos(33  333°)6378.137km  =  3.827  x  10  '  km  (162) 
Assuming  that  these  values  are  equivalent  to  the  semi-axes  of  an  error  ellipse,  multiplying 
longitude  error  by  latitude  error  and  by  n  gives  the  area  of  the  error  ellipses: 

ellipse  area^^^  =1.019  x  102  km2  (163) 

ellipse  areafoc4lpUaeeoor  =  2.574  x  10' km2  (164) 

The  focal  plane  error  effects  are  25  times  greater  than  the  effects  of  the  satellite  position 
error  The  relative  magnitudes  agree  with  the  numerical  results,  but  the  specific  values 
obtained  here  are  greater  than  those  obtained  numerically  by  a  factor  of  about  three  The 
numerically  obtained  values  were: 

launch  ellipse  mean  areattlellIIet>otlbol)eiIor  =  5.27  x  10  3  km2  (165) 

launch  ellipse  mean  areafoc,lpUneeno,  =  8.45  x  10~2  km2  (166) 

This  difference  is  expected,  since  this  analysis  does  not  account  for  stereo  viewing,  the 
least  squares  iteration  process,  or  time  effects 

This  result  shows  that  it  is  possible  to  compute  error  ellipse  areas  using  a 
simplified  qualitative  analysis  The  relative  order  of  the  magnitude  of  the  effects  is  in 
agreement  with  predictions,  and  verifies  the  results  obtained  numerically  To  do  the  same 
analysis  for  burnout  volumes  and  impact  areas  would  require  the  addition  of  time  errors, 
and  trajectory  equations  This  increases  the  difficulty  of  this  "back  of  the  envelope" 
analysis,  and  is  not  treated  here 
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VII.  CONCLUSION 


The  TALON  SHIELD  /  ALERT  state  vector  estimation  algorithm  was  correctly 
modeled  in  MATLAB™  System  errors  were  modeled  and  introduced  into  the  algorithm 
to  determine  the  effects  upon  the  accuracy  of  the  final  results  A  simplified  qualitative 
analysis  verified  the  quantitative  results. 

The  three  questions  posed  in  the  first  chapter  are  now  answered 

1 .  The  errors  in  the  system  are  broadly  categorized  as  errors  in  time,  satellite 
position,  and  LOS 

2.  The  relative  magnitudes  of  the  error  effects  listed  largest  to  least: 

A    LOS  -  using  a  zero-mean,  5  uradian  standard  deviation  normal  error 

distribution  in  focal  plane  coordinates 
B    satellite  position  -  using  a  zero-mean,  200  meter  standard  deviation  normal 

error  distribution  in  satellite  positon  coordinates 
C.  time  -  using  a  uniform  error  distribution  between  zero  and  1  millisecond 
3    The  effects  seem  to  behave  independently,  since  the  principle  of  superposition 
works    The  effects  of  the  combined  errors  are  slightly  greater  than  the  sum  of 
the  individual  effects 

Additional  runs  were  made  to  determine  the  effects  of  various  geometries,  spin  rates,  and 
observation  synchronization 

There  are  many  tasks  left  to  accomplish  in  this  area  of  research  The  MATLAB 
code  can  be  optimized,  additional  parameters  and  system  changes  can  be  simulated,  and 
more  detailed  analysis  of  the  results  can  be  made  The  error  sources  can  be  modeled  more 
exactly  by  analytically  "peeling  the  error  onion",  or  by  using  actual  DSP  values  and 
measurements    The  results  obtained  could  be  further  analyzed,  possibly  to 

Simulating  the  DSP  system  with  computer  code  has  proven  to  be  an  effective  tool 
for  error  analysis,  which  is  a  necessary  first  step  for  determining  how  best  to  improve  the 
present  TBM  detection  system  or  modify  the  design  of  a  follow-on  system 
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APPENDIX  A.  "A"  MATRIX  DERIVATION 

The  A  matrix  is  made  up  of  partial  derivatives  of  focal  plane  coordinates  with 
respect  to  the  tactical  parameters  These  elements  can  be  determined  qualitatively  Exact 
equations  are  not  required  for  the  iteration  process  to  converge,  so  small-angle 
assumptions  are  made  to  simplify  the  resulting  equations  The  benefit  of  using  these 
approximations  for  the  matrix  elements  is  the  ease  of  computation  and  thus  a  reduction  of 
computing  time 

The  A  matrix  is  divided  into  three  components,  Ai,  A2,  and  A3  The  elements  of 
Ai,  and  A;  will  be  derived  in  order,  using  MATHCAD  5  0  Plus  The  A3  matrix  is  derived 
in  Chapter  3    The  elements  of  the  first  column  of  Ai  are  partials  of  latitude,  <t> 

4-        acos  cos(G)  sin  «!►  q    -  sin(0)  cos  <f>  q    cos,  a  q   j 

Taking  the  partial  derivatives 

dp  sin(9)  sin  0  q    -  cos(0)  cos  4>  0    cos  a  0 

dG~  "  ~~2 

1       cos(9)sin0Q       sin(9)  cos  0  q    cos  a  q 

Using  the  small  angle  assumption  that  sin(0)=O  and  cos(0)=l 

dx     cos  4  0   cos  a0 

=      —I-         _  '  =cos  a  0 

do      r  2  ° 

_!l     sm>0 

e-  d 

refi 
d9      1 

s 

dd    refl 

d=(l  -  1.5L)dp 

dp«a0  +  a1  (T-T0    ^a2    T     TQ     +a3    T    TQ     +  a4    T     TQ 
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—=(1      1.5  L)  J -a,  -2a2lT-T0)-  3a3    T-T0  2-4a4    T     TQ  3 
dT0 

dd  f  2  3  4 

— -  1.5  ;a0     a  j    T    TQ  j  +  a2  J-  TQ!    +a3    T-  TQ     +a4  jT-Tq 
dL 

dp  _d+  d6    dd  _<*>*[*  0}  jid 
dTo'd0  dd  dT0"     reff      dT0 

d>_d*  d9  dd_   ]  5dpcos,a0: 
dL*  d9  dd  dL~ 


reflF 


do  0 

-dU0 
d*,o 


d+ 

—  ■  6  sin  a  n 
da0 

The  second  column  is  composed  of  partials  of  longitude,  k: 


sin(6)  sin  a  q 

/.=  /.  n  ■  asin 

cos(^) 

d'K 

1                               sin  a  o      sin  a  q 

d9     ; 

2               cos(<fr)      cos(<|>) 

jl 

7  sin  a  r. 
sin(G)2              2 
cos(o) 

'J 

d/.  _d/.  d9    dd         sm,a0         dd  _     sin,tt  0;         dd 
dT0    dG  dd  dT0    re^-cos(o)  dTQ    r^cosp  q    dT  0 

Assuming  (J>  is  approximately  <J>0 
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d/._d/.  d9  dd_   1  ^dpsina0  ^   1.5dpsina0, 
dL    d9  dd  dL        re^cos(<!))  reg-cos  *  q 

1  sin  a  Oy  do  sin(o) 

- ■  -sin(G)-       — -sii<4)  —  «an(G)  sin.  a  0  =0 

cos(o)2  d*  0  cos(o)2 


1     sin(G)' 


cos(4)2 


61 
din 


61 

dhn 


d>.  _ 
da  0  = 


cos  a  o      0  cos  a  q 
;z=  sin(9)  = 

cos(6)       cosv4  0 


1      sin(6) 


;  0 


cos(^) 
And  the  third  column  elements  are  the  partials  of  height,  h 
h=h0     (1  -  L)hp 


_dA 
dTn 


0         D3     J       l0 


■(1*L)      b,  -2-b2-T     T0       3br  T     T0        4b4    T     Tc 


dh_ 
dL~ 


-hi 


*--0 

do  n 


dh 
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— =0 
da  n 


dhf 


*0 

The  A2  matrix  is  the  partials  of  (XYZ)  with  respect  to  (<|>,  X,  h): 
x7k  =  [r^ik  cos(<J>\  )  +  ahk  cos(<|>k  )]cos(?.k ) 
yik  =  [^u  cos(4»'k  )  +  altk  cos((j>k)]sin(?.k) 
zt,  =  rioc*ik  sin(f  k  )  +  altk  sin($k) 
Approximating  r^ai  with  reff ,  alt  with  h,  $'  with  4>,  and  taking  the  partials 

dx 

—  =  -<rcff  +  h)sin(<t>)cos(),) 

dx 

—  =  -(rcff+h)cos((j))sin(X) 
cK 

—  =  cos(<|))cos(a) 
dh 

dy 

T7rr-(rcff  -h)sin(<j))sin(A) 

dy 

dy 

-f-  =  (refl  +h)cos((J))cos(X) 


o 

rh 

=  cos(<t))sin(>L) 

=  (reff  +  h)cos(4>) 

*=0 

—  -  sin((|)) 
cm 
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The  geocentric  longitude  and  latitude  can  be  expressed  in  terms  of  (XYZ) 
coordinates 


dx 


dx 


?.\\^y 


2,1 


X  •    1  +  i- 

dX  1 


2\1 


x  •   1-t 


4'=alaii 


dx 


W^+y2 
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2    ll  +  — ! 
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dz 
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The  differentials  of  the  geocentric  longitude  and  latitude  can  be  determined  using 

the  Chain  Rule    Using  longitude  and  satellite  radius  as  example  quantities 

dX  _dX    dx     dX    dy 
dR  ~  dx    dR  +  dy    dR 

Substituting  the  following  quantities  into  the  previous  equations: 

R  =42164. 17km 

gfaa   =70deg 

8  =0 

'local  =6378.1371^ 

ti  =  varies  0r  to  8.5° 
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p  =  varies  0°  to  360° 

This  gives  expressions  for  the  variations  of  the  geocentric  longitude  and  latitude 
due  to  errors  in  LOS  projections  from  a  DSP  satellite  with  a  Greenwich  hour  angle  of  70°. 
These  expressions  are  transformed  into  plots  for  examination. 

The  following  plots  show  the  differentials  of  geocentric  longitude  and  latitude  with 
respect  to  the  five  error  sources  The  horizontal  axis  pointing  to  the  left  is  rj,  the 
horizontal  axis  pointing  toward  the  right  is  0,  and  the  vertical  axis  is  the  differential 
evaluated  at  each  (r),  0)  coordinate 
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APPENDIX  C.  MATLABCODE 

The  following  is  the  MATLAB™  program  used  to  simulate  the  TALON  SHIELD 

/  ALERT  state  vector  estimation  algorithm    It  is  commented  in  blue  font 

°o  Martin  Beaulieu 

°-o  DSP  Estimation  Algorithm 

°o  Initialize  variables 


clear, 

format  long, 

numpoints=1000, 

reff=6371, 

f=l/298.257, 

re=6378.137, 

a0=0  0749291380897402, 

a  1  =-0.050 1647424642791, 

a2=0.004 194795 1949387, 

a3=-0.0000141392516462238; 

a4=8.5 1 293969732759E-07, 

b0=0. 884282320240989, 

bl=-0.129163814041924, 

b2=0.01 37288843548327, 

b3=-0.0001 5930728876899 1 , 

b4=  1 .36938557528 1 99E-06, 

tmax=62.5, 

mu=398601  2, 

d2r=pi/180, 

r2d=180/pi, 

we=15*d2r/3600, 


%  clear  memory 

%  use  1  5-digit  numerical  format 

°o  number  of  data  points  to  run 

°o  spherical  earth  radius  (km) 

°o  WGS-84  flattening  factor  for  oblate  earth 

°o  oblate  earth  equatorial  radius  (km) 

°o  coefficients  in  range  profile  equation 


i  coefficients  in  height  profile  equation 


°o  TBM  maximum  motor  burn  time  (seconds) 

°o  earth  gravitational  constant  (km  '3  sec '2) 

°o  degrees  to  radians 

°o  radians  to  degrees 

%  earth  rotation  rate  (  1  5  degrees  hour) 


%  This  section  generates  the  observation  table  (time.  az.  el.  gha.  dec.  and  radius)  The  first 
%  step  is  to  choose  the  tactical  parameters,  then  calculate  hov\  that  missile  profile  would 
°o  appear  to  an  orbiting  sensor 


°o  Initial  tactical  parameters 

for  ii=0: 1 1  °  < 

T0=100, 
L=0, 
%latO=12.783*d2r,lonO=45.050*d2r, 


loop  through  12  headings  around  compass  rose 
°o  launch  time  of  day  in  seconds 
°o  loft  parameter 
°o  Aden.  Yemen 
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Iat0=33.333*d2r,lon0=44.433*d2r, 

%latO=33.500*d2r,lonO=36.319*d2r, 

%latO=24.633*d2r,lonO=46.717*d2r, 

%latO=35.666*d2r,lonO=51.433*d2r, 

hO=0; 

hdgO=ii*pi/6, 

data=[TO  L  lonO  latO  hO  hdgO  zeros(l,9)], 


°/b  Baghdad.  Iraq 

%  Damascus,  Syria 

°-o  Riyadh.  Saudi  .Arabia 

%  Tehran.  Iran 

%  launch  altitude 

%  launch  heading  in  30  degree  steps 

%  true  tactical  parameters 


°/o  Generate  observation  time  matrix 

obsl=TO+30+5*randn(l), 

dt  1  =  1 0*rand(  1  ),dt2=  1 0*rand(  1 ), 

ttime=[obs  1  ,obs  1  +dt  1  ,obs  1  +dt2], 

ttime=sort(ttime), 

nextobs=obsl  +  10, 

while  nextobs  <=  tmax+TO, 
ttime=[ttime,ttime(j)+10], 

rh-i; 

nextobs=ttime(j)+10, 
end 
n=length(ttime), 

%  Generate  satellite  position  matrix 


>  first  observation  time  -  clouds,  vapor,  etc 

>  two  random  time  intervals.  0<dt<10  sec 

>  first  3  observation  times 

>  put  times  in  chronological  order 
)  next  sequential  observation 

>  counting  index 

)  if  next  observation  is  during  motor  burn 
)  add  next  observation  to  time  matrix 
)  increment  counter 

>  next  sequential  observation 

>  loop  until  matrix  is  completed 
,  number  of  observations 


gha=[10*d2r,70*d2r;105*d2r]; 
satord=ceil(6*rand(  1 )), 
if  satord  =  1 

tgha=[gha(l),gha(2),gha(3)], 
elseif  satord  =  2 

tgha=[gha(l),gha(3),gha(2)], 
elseif  satord  =  3 

tgha=[gha(2),gha(l),gha(3)], 
elseif  satord  =  4 

tgha=[gha(2)>gha(3),gha(l)], 
elseif  satord  =  5 

tgha=[gha(3),gha(l),gha(2)], 
elseif  satord  =  6 

tgha=[gha(3),gha(2),gha(l)], 
end 

tdec=[0,0;0], 

tradius=[42164.17;42164.17;42164.17]; 
for  i=4:n 

tgha=[tgha,tgha(i-3)], 

tdec=[tdec;tdec(i-3)], 


i  Greenwich  hour  angle 

)    randomK  select  the  order  in 

)    which  satellites  observe  TB\1 


°-o  declination 

°-o  geosynchronous  radius 

°o    make  satpos  matrix  the  same  size  i 
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tradius=[tradius,tradius(i-3,:)]; 
end 
satpos=[tradius  *cos(tdec)  *cos(tgha) 

tradius  *cos(tdec)  *sin(tgha) 

tradius  *sin(tdec)],  %  satellite  position 

%  Apply  profile  to  observation  times: 

t=ttime-T0,  %  time  of  flight 

dp=aO-t-al  *t+a2*t.  A2+a3*tA3+a4*tA4,  %  nominal  profile  distance 

d=(l-l  5*L)*dp,  %  loft-modified  distance 

hp=bO+bl  *t+b2*t.A2+b3*t  A3+b4*t  A4,  %  nominal  profile  height 

h=(l+L)*hp,  °o  loft-modified  height 

%  Generate  target  position: 

theta=d/reff,  %  earth-central  angle 
lati=pi/2-acos(cos(theta)*sin(latO)+sin(theta)*cos(latO)*cos(hdgO));°o  latitude 

Iong=lon0-t-asin(sin(theta)*sin(hdg0)./cos(lati)),  %  longitude 

alt=hO+h,  %  altitude 

gcl=atan((l-f)A2*tan(lati)),  °o  geocentric  latitude 

rloc=re*(l-f)/sqrt((l-f)A2*cos(gcl)  A2+sin(gcl).A2),  %  local  radius 
tgtpos=[(rloc  *cos(gcl)+alt  *cos(lati))  *cos(long) 

(rloc.*cos(gcl)+alt  *cos(lati))  *sin(long) 

rloc  *sin(gcl)+alt.*sin(lati)],  %  target  position 

°o  Generate  line-of-sight  \ector  and  transform  into  focal  plane  coordinates 

delta=tgtpos-satpos,  %  Ime-of-sight  vector 

r3=[],r3t=[],UEN=[],  %  reset  variables 

for  i=l:n 

rotl=[cos(tgha(i))  -sin(tgha(i))  0,sin(tgha(i))  cos(tgha(i))  0,0  0  1], 

rot2=[cos(tdec(i))  0  -sin(tdec(i));0  1  0,sin(tdec(i))  0  cos(tdec(i))], 

rot3=rotl*rot2, 

r3r=[r3t,rot3'], 

r3=[r3,rot3], 

UEN(i,:Hr3t(3*i-2:3*i,:)*delta(i,:)')';  %  Up-East-North  coordinates 

end 

tbeta=rem(-pi/2-atan2(UEN(:,3),UEN(:,2))+(2*pi),(2*pi));0o  true  azimuth 
teta=atan(sqrt(UEN(:,2)  A2+UEN(:,3)A2)./(-UEN(:,l))),    %  true  elevation 


°o  I  sing  generated  data  (ttime,  tbeta,  teta.  tgha.  tdec.  tradius), 

°o  compute  tactical  parameters  using  least-squares  iteration  method 
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for  index=l  numpoints, 
index=round(index), 


%  ensure  index  is  an  integer 


o  Initialize  variables 

clear  time  radius  gha  dec  beta  eta  TO  L  latO  lonO  hO  altO  hdgO  tp, 

sse=l,dsse=l,ssel=l,failindex=0,r3t=[],r3=[],  %  reset  variables 

errtime=zeros(n,l);  %  time  error 

errrad=zeros(n,l),  %  satellite  radius  error 

errgha=zeros(n,l),  %  satellite  gha  error 

errdec=zeros(n,l),  %  satellite  dee  error 

errbeta=zeros(n,  1 ),  %  azimuth  error 

erreta=zeros(n,  1 ),  %  elevation  error 


)  Change  time,  beta.  eta.  and  satpos  by  some  error,  epsilon 


if  index  >  1 

errtime=le-3*(rand(n,l)-0  5), 

errrad=0.2*randn(3,l), 

errgha=atan(0.2/42 1 64  1 7)*randn(3, 1 ), 

errdec=atan(0.2/42 1 64  1 7)*randn(3, 1 ); 

for  i=4:n 

errTad=[enTad,errrad(i-3)], 
errgha=[errgha,errgha(i-3)], 
errdec=[errdec,errdec(i-3)], 

end 

errbeta=5e-6*randn(n,  1 ); 

erreta=5e-6*randn(n,  1 ), 
end 

time=ttime+errtime, 
radius=tradius+errrad, 
gha=tgha+errgha, 
dec=tdec+errdec, 
beta=tbeta+errbeta, 
eta=teta+erreta, 


o-o  <=--0.00]  second 
°o  variance  of -'--200  meters 
°o  variance  of -'--200  meters 
°o  variance  of --200  meters 
°o  fix  satellite  positions 
°o        during  sinule  run 


°o  variance^  -omicroradians 
°o  variance=  -omicroradians 

°o  time  =  truth  -  error 
°o  radius  =  truth  -  error 
°o  gha  =  truth  +  erroi 
°o  dec  =  truth  -  error 
°o  beta  =  truth  -  error 
°o  eta  =  truth  ~  error 


°-o  Compute  target  position 

satpos=[radius  *cos(dec).  *cos(gha) 
radius.  *cos(dec).*sin(gha) 

radius. *sin(dec)],  °o  satellite  position 

for  i=l:n  %  rotation  matrices 

rotl=[cos(gha(i))  -sin(gha(i))  0;sin(gha(i))  cos(gha(i))  0,0  0  1]; 
rot2=[cos(dec(i))  0  -sin(dec(i));0  1  0,sin(dec(i))  0  cos<dec(i))], 
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rot3=rotl*rot2, 
r3t=[r3t,rot3*], 
r3=[r3,rot3], 
uen=[-cos(eta(i)), 

-sin(beta(i))*sin(eta(i)), 

-cos(beta(i))*sin(eta(i))], 
ehat=rot3*uen, 

chi=pi-asin((norm(satpos(i,:))*sin(eta)/reflr)); 
rho=reff7sin(eta)  *  sin(eta+chi), 
los=rho*ehat, 
tgtpos=los+satpos(i, )', 

lat(i)=atan(tgtpos(3)/sqrt(tgtpos(  1  )A2+tgtpos(2)A2)): 
lon(i)=atan2(tgtpos(2),tgtpos(  1 )), 
end 


°/o  Up-East  -Nonh  coordinates 

%  line  of  sight  direction 

°o  computation  of  length 

°o  for  los  vector 

°-o  line-of-sight  vector 

0  o  target  position 

o  latitude 

o  loncitudc 


°o  Intial  estimate  of  tactical  parameters 


latO={lat(l)+lat(2)+lat(3))/3, 
lonO=(lon(l)+lon(2)+lon(3))/3, 
latbo=(lat(n-2)+lat(n-l)+lat(n))/3, 
lonbo=(lon(n-2)+lon(n- 1  )+lon(n))/3 , 


°o  first  ohs  latitude 
°o  first  obs  longitude 
°o  latitude  at  last  obs 
%  lonmtude  at  last  obs 


hdgO=rem(pi/2-atan2(latbo-latO,(lonbo-lonO)*cos(latO))+(2*pi),(2*pi)),  °  olneh  heading 

T0=time(l)-20,  %  launch  time 

L=0,  °o  loft  parameter 

h0=0,  °  o  launch  height  abo\  e  WGS-84  ellipsoid 

tp=[TO,L,latO;JonO,hO,hdgO],  °o  tactical  parameters  matrix 


%  Iterate  on  initial  estimates  until  the  difference  between  the  sum  of  the 
%  squares  of  the  errors  between  two  consecutive  runs  is  <=10*eps: 


while  dsse  >  10*eps 
duv=[],A=[], 

TO=tp(l),L=tp(2),latO=tp(3), 
Ion0=tp(4),h0=tp(5),hdg0==rem(tp(6),2*pi), 
t=time-T0; 

dp=aO+al*t+a2*t.A2+a3*t.A3+a4*t.A4, 
d=(l-l  5*L)*dp, 

hp=bO+bl*t+b2*tA2+b3*t.A3+b4*tA4, 
h=(l+L)*hp, 
theta=d/reff, 

lati=pi/2-acos(cos(theta)*sin(latO)+sin(theta)*cos(latO)*cos(hdgO)),0  b  geodtic  latitude 
long=lonCH-asin(sin(theta)*sin(hdgO)/cos(lati)),  °o  longitude 
alt=hO+h,  °o  altitude 

gcl=atan((l-f)A2*tan(lati)),  %  geocentric  latitude 

rloc=re*(l-f)/sqrt((l-f)A2*cos(gcl).A2+sin(gcl)  A2),°o  local  radius 


°o  loop  until  dsse  <  10*eps 

°o  reset  \ariables 

°o  update  tp  matrix  values 

°o  with  dtp's  added 

0  o  time-of-flight 

°o  nominal  profile  distance 

%  loft-modified  distance 

%  nominal  profile  height 

°o  loft-modified  height 

°o  earth-central  ande 
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tgtpos=[(rloc  *cos(gcl)+alt  *cos(lati)).*cos(long) 
(rloc. *cos(gcl)+ah  *cos(lati)). *sin(long) 

rloc.*sin(gcl)+alt  *sin(lati)];  %  target  position 

delta=tgtpos-satpos,  °o  line-of-sight  vector 

for  i=l  :n  %  rotate  into  UEN  coordinates 

UEN(i,:)=(r3t(3*i-2:3*i,:)*delta(i,:)')', 
end 

uhat=-UEN(:,2)./UEN(:,l);  %    estimated  focal 

vhat=-UEN(:,3)./UEN(:,l),  %    plane  coordinates 

u=-tan(eta).*sin(beta),  °o         actual  focal 

v=-tan(eta)  *cos(beta),  %         plane  coordinates 

du=u-uhat,  °o  difference  between 

dv=v-vhat,  %  estimate  and  actual 

°o  Compute  A  matrix: 

cl=cos(hdgO)/reff,c2=sin(hdgO)/reff7cos(latO),     °o  constants  in  A  matrix 
c3=sin(hdgO)/reff,c4=cos(hdgO)/reff7cos(latO), 

dddT0=-(l-1.5*L)*(al+2*a2*t+3*a3*t.A2+4*a4*tA3),  %  -horizontal  speed 
dhdT0=-(l+L)*(bl+2*b2*t+3*b3*tA2+4*b4*tA3),        %  -vertical  speed 
dudUEN=[UEN(:,2)./UEN(:,l)A2  -l./UEN(:,l)  zeros(size(time))], 
dvdUEN=[UEN(:,3)./UEN(:,l)  A2  zeros(size(time))  -l./UEN(:,l)]; 
°o  change  in  focal  plane  coordinates  w  rt  changes  in  I  EN  coordinates 
for  i=l:n  °o  construct  A  matrix 

duv=[duv,du(i),dv(i)],  "  o  error  \  alues  matrix 

Al=  [dddT0(i)*cl  dddT0(i)*c2  dhdTO(i), 
-1.5*dp(i)*cl  -1.5*dp(i)*c2  hp(i), 
1  0  0,0  1  0,0  0  l,-d(i)*c3  d(i)*c4  0]; 
A2=[-(rloc(i)+alt(i))*sin(lati(i))*cos(long(i)) 
-(rloc(i)+alt(i))*sin(lati(i))*sin(long(i)) 
(rloc(i)+alt(i))*cos(lati(i)), 
-(rloc(i)+alt(i))*cos(lati(i))*sin(long(i)) 
(rloc(i)+alt(i))*cos(lati(i))*cos(long(i))0, 
cos(lati(i))*cos(long(i))  cos(lati(i))*sin(long(i))  sin(lati(i))], 
Atmp=Al*A2*r3(3*i-2:3*i,:),  %  A3=(delta  UEN  /  delta  \yz)=r3 

A=[A;(Atmp*dudUEN(i,:)7,(Atmp*dvdUEN(iv )')'],  °  o  total  A  matrix 

end 

AT=pinv(A),  °  o  find  pseudoinverse  of  A 

dtp=AT*duv,  °o  find  changes  to  tp 

tp=tp+dtp;  °o  add  changes  to  tp 

ssel=sum(duvA2),  %  find  sum  of  the  squares  of  the  error 

dsse=abs(sse-ssel);  %  difference  between  iterations 

if  ssel  <  sse 
sse=ssel; 
end 
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failindex=failindex+ 1 ,  °o  increment  counter 

if  failindex>  1 00  %  if  statement  to  prevent  continuous  looping 

dsse=0,  °o  in  the  event  that  it  does  not  converge 
end 

end  %  loop  to  next  tp  iteration 

°o  Burnout  time  estimation 

Tlast=time(n),  %  last  observation  time 

Tnext=2*time(n-2)-time(n-5),  %  next  observation  if  TBM  was  still  burning 

Tmax=T0+tmax,  %  max  observation  time  from  profile 
if  index  =1 

Tbo=T0+tmax, 

elseif  Tmax>=Tnext  %  pick  burn  out  time  as  half  of  the  time 

Tbo=Tlast+0.5*(Tnext-Tlast),  %  between  Tlast  and  Tnext  or  Tlast  and 

else  °o  Tmax  depending  on  relative 

Tbo=Tlast-K).5*(Tmax-Tlast),  %  magnitude  of  Tmax  and  Tnext 
end 

tbo=Tbo-T0;  %  burnout  time-of-flight  (close  to  tmax) 

°o  Calculate^  state  \ector  at  burnout 

d=(l-1.5*L)*(a(Hal*tbo+a2*tboA2+a3*tboA3+a4*tboA4),  °o  distance 

ddot=(l-l  5*L)*(al+2*a2*tbo+3*a3*tboA2+4*a4*tboA3),  %  horizontal  speed 

h=(l+L)*(b0+bl*tbo+b2*tboA2+b3*tboA3+b4*tboA4),  %  height 

hdot=(l+L)*(bl+2*b2*tbo+3*b3*tboA2+4*b4*tboA3),  °/o  vertical  speed 

theta=d/reff,  %  earth-central  angle 
latbo=pi/2-acos(cos(theta)*sin(latO)+sin(theta)*cos(latO)*cos(hdgO));0ogeodtic  latitude 

Ionbo=lonO+asin(sin(theta)*sin(hdgO)/cos(latbo)),  %  burnout  longitude 

hbo=h0+h,  °o  burnout  altitude 

Vbo=sqrt(ddotA2+hdotA2),  %  burnout  velocity 

fpabo=atan(hdot/ddot),  %  flight  path  angle 

hdgbo=rem(asin(cos(latO)*sin(hdgO)/cos(latbo))+(2*pi),(2*pi)),  %  burnout  heading 

°o  Calculates  ballistic  trajectory  of  target  to  impact 

altbo=h0+hbo,  %  burnout  altitude 

gclbo=atan(((  1  -f)A2)*tan(latbo)),  °  o  geocentric  latitude 

rlocbo=re*(  1  -f)/sqrt((  1  -f)A2*cos(gclbo)A2+sin(gclbo)A2),  %  local  radius 

tgtposbo=[(rlocbo*cos(gclbo)+altbo*cos(latbo))*cos(lonbo) 
(rlocbo  *cos(gclbo)-t-altbo  *  cos(latbo))*  sin(lonbo) 
rlocbo*sin(gclbo)+altbo*sin(latbo)];    %  TBM  position  at  burnout 
rbo=norm(tgtposbo),  °0  radius  at  burnout 

vo=rlocbo*we*cos(gclbo),  %  initial  inertia!  velocitv 

Vu=Vbo*sin(fpabo),  %  up  component  of  velocit\ 
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%  east  component  of  velocity 

°o  north  component  of  velocity 

%  inenial  speed 

%  inenial  flight  path  angle 

%  inertial  heading 

%  energy  parameter 

%  eccentricity  of  orbit 

%  freeflight  eca 

%  eccentric  anomalv 


Ve= Vbo  *  cos(fpabo)  *  sin{hdgbo)+vo , 
Vn= Vbo  *  cos(fpabo)  *  cos(hdgbo), 
Vin=sqrt(Vu.A2+Ve.A2+Vn.A2), 
fpain=asin(Vu/Vin), 

hdgin=rem(pi/2-atan2(-Ve,Vin)+2*pi,2*pi), 
Qin=VinA2*rbo/mu, 
ein=sqrt(  1  +Qin*(Qin-2)*cos(fpain)A2); 
PSI=2*acos((  1  -Qin*cos(fpain)A2)/ein)+theta, 
E=acos((ein-cos(PSI/2))/(  1  -ein*cos(PSI/2))), 
ain=rbo/(2-Qin), 

tflF=2*sqrt(ainA3/mu)*(pi-E+ein*sin(E)), 
gclirn=asin(siri(gclbo)*cos(PSI)+cos(gclbo)*sin(PSI)*cos(hdgin)),  %geocentric  latitude 
latim=atan(tan(gclim)/(l-f)A2),  %  geodetic  latitude  of  impact 

fiHon=acos((cos(PSI)-sm(gclim)*siji(gclbo))/(cos(gclim)*cos(gclbo)))-we*tff, 
lonim=lonbo+fflon,  %  longitude  of  impact 

tim=2*tbo+tff,  °-o  total  flight  time 

Tim=TO+tim,  %  time  of  day  of  impact 

data=[data;tp(l)  tp(2)  tp(4)  tp(3)  tp(5)  tp(6)  lonbo  latbo  hbo  fpabo  hdgbo  Vbo  Tim 
lonim  latim],  %  sa\  e  data  in  matrix 

end  °o  end  of  numpoint  loop 


%  semi  major  axis 
%  time  of  free  fliaht 


°o  Plottinu  routines  and  data  analvsis 


tllon=data(2,3)*r2d; 
tllat=data(2,4)*r2d, 
Uon=data(3 :  index+ 1 ,3  )*r2d, 
llat=data(3 :  index+ 1 ,4)*r2d, 
mllon=mean(data(3  :index+ 1 ,3))*r2d, 
mllat=mean(data(3 :  index+ 1 ,4))*r2d, 
dllon=tllon-mllon, 
dllat=tllat-mllat, 


°o  true  launch  longitude 

°o  true  launch  latitude 

°o  corrupt  launch  longitude 

%  cornipt  launch  latitude 

°o  mean  of  corrupt  longitudes 

°-o  mean  of  corrupt  latitudes 

°o  difT between  mean  and  true  lion 

°o  diff  between  mean  and  true  Hat 


tblon=data(2,7)*r2d, 
tblat=data(2,8)*r2d; 
tbalt=data(2,9), 
blon=data(3 :  index+ 1 ,7)*r2d, 
blat=data(3  :index+ 1 ,8)*r2d, 
balt=data(3 :  index+ 1 ,9), 
mblon=mean(data(3 :  index+ 1 ,7))*r2d, 
mblat=mean(data(3 :  index+ 1 ,8))*r2d, 
mbalt=mean(data(3 :  index+ 1 ,9)), 
dblon=tblon-mblon, 
dblat=tblat-mblat, 
dbalt=tbalt-mbalt, 


j  true  burnout  longitude 

)  true  burnout  latitude 

a  true  burnout  altitude 

d  corrupt  burnout  longitude 

d  conupt  burnout  latitude 

d  corrupt  burnout  altitude 

3  mean  of  corrupt  bo  longitudes 

3  mean  of  corrupt  bo  latitudes 

3  mean  of  corrupt  bo  altitudes 

3  difT  between  mean  and  true  blon 

o  ditT between  mean  and  true  blat 

o  difT  between  mean  and  true  blat 
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tilon=data(2,14)*r2d, 
tilat=data(2,15)*r2d, 
ilon=data(3 :  index+ 1 , 1 4)*r2d, 
ilat=data(3  :index+ 1 , 1 5)*r2d, 
milon=mean(data(3  :index+ 1 , 1 4))*r2d, 
milat=mean(data(3 :  index+ 1 , 1 5))*r2d, 
dilon=tilon-milon, 
dilat=tilat-milat, 


°-o  true  impact  longitude 

%  true  impact  latitude 

%  corrupt  impact  longitude 

°o  corrupt  impact  latitude 

%  mean  of  corrupt  im  longitudes 

°o  mean  of  corrupt  im  latitudes 

°o  diff  between  mean  and  true  lion 

°'o  diff  between  mean  and  true  ilat 


lpts=pion-tllon  llat-tllat], 

lk=cov(lpts), 

[lv,ld]=eig(lk), 

lab=2*sqrt(diag(ld)), 

la=lab(l), 

lb=lab(2), 

lth=pi/2-atan(lv(  1 , 1  )/lv(2, 1 )), 


°.o  center  data  about  launch  truth 

%  covariance  of  scatter 

°o  eigen  value  &.  vector  of  covariance 

°o  dimensions  of  ellipse 

°o  launch  ellipse  semimajor  axis 

%  launch  ellipse  semiminor  axis 

°o  angle  to  rotate  ellipse 


bpts=[blon-tblon  blat-tblat  balt-tbalt], 
[bvl,bdl  ]=eig(cov(bpts(,  1  ),bpts(:,2))), 
babc=2*sqrt(diag(bd  1 )), 
ba=babc(  1  ),bb=babc(2),bc=babc(3), 


°o  center  data  about  burnout  truth 
°o  eigen  val  and  vector  of  co\ariance 
°o  dimensions  of  ellipse 
°o  burnout  ellipse  semimaior  axis 


ipts=[ilon-tilon  ilat-tilat]; 

ik=cov(ipts), 

[iv,id]=eig(ik), 

iab=2*sqrt(diag(id)), 

ia=iab(l), 

ib=iab(2), 

ith=pi/2-atan(iv(  1 , 1  )/iv(2, 1 )), 


°o  center  data  about  launch  truth 

°c  co\  ariance  of  scatter 

°o  eigen  \alue  k  \ector  of  coxarance 

°(.  dimensions  of  ellipse 

°o  impact  ellipse  semimajor  axis 

°o  impact  ellipse  semiminor  axis 

°o  angle  to  rotate  ellipse 


°o  Sa\e  data  to  anahze 

data=[data,la  lb  1th  bal  bbl  bthl  ba2  bb2  bth2  ba3  bb3  bth3  ia  ib  ith,  mllon  mllat  mblon 
mblat  mbalt  milon  milat  dllon  dllat  dblon  dblat  dbalt  dilon  dilat  0], 
ifu==0 

save  bOt  data  -ascii  -double 
elseif  ii=l 

save  b3t  data  -ascii  -double 
elseif  ii=2 

save  b6t  data  -ascii  -double 
elseif  ii=3 

save  b9t  data  -ascii  -double 
elseif  ii=4 

save  bl2t  data  -ascii  -double 
elseif  ii=5 
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save  bl  5t  data  -ascii  -double 
elseif  ii=6 

save  b  1 8t  data  -ascii  -double 
elseif  ii=7 

save  b21t  data  -ascii  -double 
elseif  ii=8 

save  b24t  data  -ascii  -double 
elseif  ii=9 

save  bt27  data  -ascii  -double 
elseif  ii= 10 

save  bt30  data  -ascii  -double 
elseif  ii=l  1 

save  bt33  data  -ascii  -double 
end 

clear  data 
end 
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